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Measurement  and  Calculation  of  Developing  Turbulent  Flow 
in  a  U-Bend  and  Downstream  Tangent  of  Square  Cross-Section 


ABSTRACT 

Experimental  measurement  and  numerical  modeling  has  been  performed 
on  the  flow  in  a  180°  bend  of  square  cross-section,  preceeded  and  follow¬ 
ed  by  straight  ducts  of  the  same  cross-section. 

Measurements  of  mean  velocities  and  their  associated  turbulent 
stresses  along  the  streamwise  (e)  and  the  gapwise  or  radial  (r)  direc¬ 
tions  have  been  made  at  several  longitudinal  planes,  using  the  non- 
intrusive  laser-Doppler  velocimeter  technique  in  backscatter  mode. 

Mean  flow  data  reveal  features  in  qualitative  agreement  with  results 
obtained  from  inviscid  flow  analysis.  Measurements  of  the  turbulent 
stresses  display  previously  undocumented  anisotropic  characteristics 
arising  from  shearing  motions  induced  in  the  core  of  the  flow.  In  the 
downstream  tangent,  measurements  show  that  drastic  reductions  of  the 
secondary  motion  take  place  in  less  than  5  hydraulic  diameters.  From 
then  on,  however,  the  flow  recovers  only  slowly  from  the  effects  of 
the  bend. 


Numerical  predictions  using  a  k-e  model  of  turbulence  have  been 
performed  for  the  experimental  configuration.  Model  deficiencies  are 


more  clearly  revealed  by  the  higher  order  accuracy  of  the  finite  differ¬ 
encing  scheme  used  and  the  implementation  of  a  partially-parabolic  cal¬ 
culation  algorithm.  It  is  found  that,  with  the  problem  of  numerical 
diffusion  relieved,  false  physical  diffusion  and  the  isotropic  character¬ 
istics  inherent  in  the  model  are  the  main  causes  for  the  differences 
observed  between  measurements  and  computations. 

Additional  predictions  with  an  algebraic  stress  model  (ASM)  show 
poor  agreement  with  the  measurements.  As  a  result,  two  "experimental" 
tests  have  been  carried  out  to  check  indirectly  the  predicting  ability 
of  the  ASM  closure.  It  is  found  that  the  proposed  model  is  capable  of 
resolving  the  anisotropy  of  the  turbulence  correctly,  provided  that  the 
mean  velocity  field  is  known  accurately. 


FOREWORD 


This  report  represents  the  culmination  of  a  four  year  effort  under¬ 
taken  at  the  Berkeley  Campus,  aimed  at  measuring  and  making  predictable 
the  turbulent  flow  in  passage  through  a  180°  curved  duct  of  square  cross 
section.  The  research  program  initially  started  out  as  a  joint  venture 
with  the  Davis  campus  of  the  University  of  California.  After  year  one, 
the  Davis  portion  of  the  program,  under  the  direction  of  Professor  Brian 
E.  Launder,  was  moved  to  the  University  of  Manchester  Institute  of  Science 
and  Technology.  In  spite  of  the  distances  and  difficulties  involved,  the 
close  and  beneficial  collaboration  with  Professor  Launder  and  his  group 
at  UMIST  continued  uninterruptedly  throughout  the  duration  of  the 
program  in  Berkeley.  We  are  indebted  to  Mr.  K.  Ellingsworth  of  the 
Office  of  Naval  Research  for  making  the  collaboration  possible. 

Whereas  research  in  UMIST  has  centered  on  the  heat  transfer  aspects 
of  turbulent  flow  in  bends,  that  in  Berkeley  has  dealt  exclusively  with 
the  fluid  mechanics.  At  the  time  of  writing  measurements  are  being 
obtained,  in  Berkeley,  in  a  180°  curved  pipe  to  extend  the  data  base 
available.  The  results  will  be  reported  in  a  subsequent  addendum  to  this 
report.  The  data  will  be  used  by  the  UMIST  team  to  test  and  validate  a 
numerical  procedure  for  predicting  turbulent  flow  in  curved  ducts  of 
circular  cross-section. 
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#1,  ( . )  "experimental"  test  #2. 


CHAPTER  1.  INTRODUCTION 


1 .1  The  Problem  of  Interest 

Considerable  efforts  have  been  devoted  to  the  study  of  flow  in 
curved  conduits  over  the  past  40  years  due  to  its  fundamental  as  well 
as  industrial  importance.  When  a  fluid  flows  through  a  curved  duct  of 
square  cross-section,  the  radial  pressure  gradient  established  in  the 
core  region  (directed  to  the  center  of  curvature)  acts  upon  the  slower 
moving  fluid  near  the  side  walls,  causing  the  fluid  there  to  move  to 
the  inner  (convex)  side  of  the  bend.  In  turn,  faster  moving  fluid  in 
the  core  region  moves  to  the  outer  (concave)  side  of  the  bend,  producing 
the  secondary  motion  depicted  in  Figure  1.1.  In  addition  to  the  pres¬ 
sure-gradient  driven  secondary  flow  described  above,  gradients  of  the 
turbulent  stresses  in  a  plane  normal  to  the  main  flow  direction  can 
also  produce  secondary  motion  in  ducts  of  non-circular  cross-section. 

The  sense  of  this  secondary  motion  in  a  straight  duct  of  square  cross- 
section  is  shown  in  Figure  1.2  over  a  quarter  of  the  cross-section.  It 
carries  high  momentum  fluid  from  the  core  region  toward  a  corner  along 
the  corner  bisector  and  pushes  low  momentum  fluid  near  the  walls  back 
to  the  central  region  along  the  wall  bisectors.  This  secondary  motion 
causes  the  contours  of  mean  longitudinal  velocity  to  bulge  toward  the 
corners  and,  as  a  result,  enhances  the  heat  transfer  rate  and  shear 
stress  there.  In  the  literature,  the  secondary  motion  induced  by  stream¬ 
line  curvature  is  usually  referred  to  as  the  "secondary  motion  of  the 
first  kind"  whereas  that  produced  by  gradients  of  turbulent  stresses  is 


sometimes  called  the  "secondary  motion  of  the  second  kind."  The  theory 
of  both  types  of  secondary  motions  has  been  reviewed  by  Johnston  [1978]. 

The  flow  configuration  of  interest  to  the  present  study  consists 
of  a  180°  circular  bend  of  square  cross-section  preceeded  and  followed 
by  straight  ducts  or  tangents.  This  U-shaped  geometry  is  found  in 
numerous  piping  systems  and  industrial  flow  configurations.  In  par¬ 
ticular,  it  is  directly  relevant  to  compact  heat  exchangers.  The  flow 
patterns  established  in  and  downstream  of  a  U-bend  strongly  influence 
the  heat  transfer  characteristics  which,  in  turn,  dictate  the  overall 
performance  of  a  heat  exchanger.  The  study  of  this  type  of  flow  also 
provides  further  understanding  of  turbulent  flow  phenomena  in  turbo¬ 
machinery  components  and,  thus,  can  lead  to  the  improved  design  of 
these  components. 

1.2  Objectives 

The  primary  objectives  of  this  study  are: 

1)  To  obtain  measurements  of  the  mean  flow  and  turbulence  charac¬ 
teristics  by  means  of  laser-Ooppler  velocimetry  in  a  U-bend  and  down¬ 
stream  tangent  of  equal  square  cross-section. 

Because  the  present  configuration  involves  the  combined  effects  of 
flow  in  a  180°  bend  and  the  downstream  tangent,  a  complex  flow  of  very 
distinct  and  previously  undocumented  velocity  characteristics  emerges. 
Besides  serving  to  advance  understanding  of  the  flow,  the  experimental 
data  should  also  be  of  value  for  developing  and  testing  turbulence 
models  suitable  for  three-dimensional  complex  turbulent  flows. 


2)  To  test  the  respective  capabilities  of  a  two  equation  (k-e) 
model  and  an  algebraic  stress  model  of  turbulence  for  predicting  the 
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mean  flow  and  its  turbulence  characteristics. 

Previous  elliptic  calculations  of  Humphrey  et  al.  [1981]  with  a 
k-e  model  for  turbulent  flow  in  a  90°  bend  of  square  cross-section 
yielded  results  in  poor  agreement  with  their  measurements  for  bend 
angles  greater  than  45°.  However,  in  that  study  it  was  not  possible 
to  separate  satisfactorily  the  inadequacy  of  the  k-e  model  from  inaccur¬ 
acies  of  numerical  origin.  Subsequent  work  by  Chang  et  al .  [1982]  using 
a  parti  ally-parabolic  procedure  and  a  quadratic  upwind -weigh ted  inter¬ 
polation  for  convection  terms  (the  QUICK  scheme)  yielded  substantially 
improved  predictions.  In  the  present  study  their  numerical  procedure 
is  employed  whenever  possible  to  establish  more  precisely  the  performance 
and  limitations  of  the  two  turbulence  models  explored. 

1.3  Earlier  Work 

1.3.1  Experimental  Studies 

Measurements  of  turbulent  flows  in  bends  of  constant  rectangular 
cross-sections  have  been  made  by  Joy  [1950],  Eichenberger  [1952,  1953] 
and  Squire  [1954]  and,  more  recently,  by  Bruun  [1979],  Humphrey  et  al. 
[1981]  and  Taylor  et  al.  [1982].  Joy  [1950]  measured  the  total  pressure 
distributions  in  three  rectangular  bends  (aspect  ratio  2:1)  of  different 
radii  to  the  bend  center  for  two  turning  angles  (either  90°  or  180°). 

He  discovered  that  for  the  180°  configuration,  the  direction  of  secondary 
motion  was  reversed  at  a  bend  angle  around  e  *  135°.  Further  experimental 


evidence  of  the  oscillatory  nature  of  secondary  motions  in  bends  of  pro¬ 
longed  curvature  were  also  found  by  Eichenberger  [1952]  and  Squire  [1954], 
and  was  summarized  by  Hawthorne  [1963]. 

In  a  study  on  the  influence  of  inlet  boundary  layer  thickness  on 
secondary  flow  development  and  pressure  losses,  Bruun  [1979]  performed 
detailed  measureements  of  mean  values  of  the  total  pressure,  static 
pressure  and  velocity  fields  as  well  as  limited  measurements  of  tur¬ 
bulence  intensity  distributions  in  two  120°  bends  of  rectangular  cross- 
sections.  By  analyzing  the  measurements,  he  was  able  to  present  a 
physical  description  of  the  secondary  flow  processes. 

The  measurements  of  developing  mean  flow  and  turbulence  charac¬ 
teristics  in  a  90°  bend  of  square  cross-section  by  Humphrey  et  al. 

[1981]  and  Taylor  et  al.  [1982]  are  directly  related  to  the  present 
work.  These  studies  were  performed  by  using  Laser-Doppler  velocimetry 
techniques  in  the  same  bend  for  the  same  Reynolds  number.  The  only 
difference  between  the  two  experiments  was  the  boundary  layer  thickness 
of  the  flow  at  the  bend  inlet  plane.  In  the  work  of  Taylor  et  al. 

[1982]  thin  boundary  layers  were  purposefully  induced  at  the  bend  inlet 
plane,  whereas  in  Humphrey  et  al.  [1981]  essentially  fully  developed 
straight  duct  flow  was  presented  as  the  inlet  condition.  The  most 
marked  differences  between  the  results  of  these  two  studies  were,  in 
general,  the:  a)  stronger  secondary  motion,  and  b)  higher  levels  of 
turbulence,  experienced  by  the  flow  with  thicker  inlet  boundary  layers. 

Detailed  measurements  of  the  much  weaker  Reynolds  stress  driven 
secondary  motions  for  the  case  of  developing  flow  in  a  duct  of  square 
cross-section  have  been  reported  by  Melling  and  Whitelaw  [1976],  and  in 


a  series  of  papers  by  Gessner  et  al.,  the  latest  reference  being  [1982]. 
It  is  noteworthy  that  the  flow  field  measurements  of  Humphrey  et  al . 
[1981]  and  Taylor  et  al.  [1982]  do  not  show  the  symmetric  corner  vortex 
structures  characteristic  of  straight  duct  flow.  This  supports  the 
notion  that  in  the  bend  sections  of  these  studies  the  cross-stream  flow 
was  influenced  more  strongly  by  the  transverse  mean  pressure-gradient 
than  by  gradients  of  the  Reynolds  stresses. 

1.3.2  Numerical  Calcul ation/Modell ing 

Curved  duct  turbulent  flow  calculations  based  on  the  numerical 
solution  of  the  Navier-Stokes  equations  have  been  made  by  Humphrey, 
Whitelaw  and  Yee  [1981]  for  their  experimental  configuration,  and  more 
recently  by  Buggeln,  Briley  and  McDonald  [1980]  for  the  flow  measured 
by  Taylor,  Whitelaw  and  Yianneskis  [1982].  In  the  work  of  Humphrey  et 
al.  [1981]  a  two-equation  k-e  model  of  turbulence  was  used  whereas 
Buggeln  et  al.  [1980]  employed  the  one-equation  model  investigated  by 
Shamroth  and  Gibeling  [1979].  In  general,  both  studies  show  good  agree¬ 
ment  between  measurements  and  calculations  of  the  mean  velocity  compon¬ 
ents  up  to  a  bend  angle  of  approximately  45°,  but  as  of  this  point  sig¬ 
nificant  deviations  appear.  For  bend  angles  larger  than  45°  calculations 
of  Buggeln  et  al.  [1980]  are  in  better  accord  with  the  measurements  than 
corresponding  results  of  Humphrey  et  al.  [1981].  This  is  attributed 
primarily  to  the  thinner  boundary  layer  inlet  condition  present  in  the 
measurements  of  Taylor  et  al.  [1982],  which  was  responsible  for  the  pro¬ 
duction  of  weaker  secondary  motions  than  were  observed  and  predicted  by 
Humphrey  et  al.  [1981].  It  should  be  remarked  that  cost  considerations 
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dictated  the  use  of  coarse  grids  in  both  of  the  above  studies  and,  since 
significant  levels  of  numerical  diffusion  probably  affected  the  results, 
there  is  an  inadequate  basis  for  judging  the  merits  and  demerits  of  the 
turbulence  models  used  respectively  by  the  authors. 

A  review  of  computational  methods  for  internal  flows  has  been  car¬ 
ried  out  by  McNally  and  Sockol  [1981].  The  following  is  a  brief  summary 
relating  to  numerical  techniques  of  relevance  to  this  work.  Within  the 
class  of  numerical  methods  which  neglects  streamwise  diffusion  is  the 
partially-parabolic  procedure  proposed  by  Pratap  and  Spalding  [1975]  for 
calculating  flow  in  curved  ducts.  In  this  method,  iterated  forward 
marching  sweeps  of  the  three-dimensional  flow  field  are  performed, 
yielding  solutions  of  the  momentum  equations  and  corrections  to  the 
pressure  field,  until  local  and  global  continuity  are  attained.  This 
type  of  numerical  procedure  was  used  by  Chang  et  al.  [1982]  to  predict 
the  flow  of  Humphrey  et  al .  [1981]  after  modifying  the  procedure  to  in¬ 
clude  quadratic  upstream  interpolation  of  convection  terms  as  in  Han 
et  al.  [1981]  to  reduce  numerical  diffusion.  In  general,  better  agree¬ 
ment  was  found  in  Chang  et  al.  [1982]  than  in  Humphrey  et  al.  [1981] 
between  mean  velocity  measurements  and  calculations  as  a  consequence  of 
the  numerical  improvements  to  the  calculation  scheme. 

Inviscid  flow  approximations  for  the  calculation  of  flows  in  curved 
ducts  have  been  proposed  by,  for  example,  Briley  and  McDonald  [1979]. 

The  numerical  procedure  developed  by  these  authors  represents  an  exten¬ 
sion  of  viscous  forward  marching  methods  in  that  it  accounts  approximate¬ 
ly  for  transverse  variations  in  streamwise  pressure  gradient.  While  this 
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approach  appears  particularly  promising  for  fast  and  relatively  inexpen- 
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sive  curved  flow  computations,  work  continues  on  its  development  and 
application;  see  Krekovsky,  Briley  and  McDonald  [1980]. 

In  contrast  to  the  pressure-dominated  flows  arising  in  curved  ducts, 
the  difficulties  associated  with  predicting  Reynolds  stress  driven  secon¬ 
dary  motions  in  ducts  of  non-circular  cross-section  are  principally  re¬ 
lated  to  the  accurate  modeling  of  cross-stream  turbulent  flow  anisotropy. 
Brundrett  and  Baines  [1964]  have  shown  that  the  turbulence  term  contribut 

ing  most  strongly  to  mean  streamwise  vorticity  in  steady,  incompressible, 
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constant  property  flow  is  (uy  "  u2  )•  In  a  later  study,  Perkins 
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[1970]  also  demonstrated  the  importance  of  the  term  {•*-*  -  -2-^)  u  u  . 

3z  3y  y  z 

In  these  terms,  y,  z,  u^  and  uz  are  the  cross-stream  coordinates  and 
velocity  fluctuations  respectively.  It  can  readily  be  shown  that  tur¬ 
bulence  models  based  on  the  notion  of  an  isotropic  turbulent  viscosity, 
such  as  the  k-e  model ,  cannot  account  for  the  generation  of  stress- 
driven  secondary  motions.  In  their  recent  study,  Naot  and  Rodi  [1981] 
have  reviewed  briefly  various  investigations  employing  Reynolds  stress 
level  closures  to  predict  these  motions.  They  conclude  that  the  essen¬ 
tial  basic  flow  features  can  be  predicted  qualitatively  with  relatively 
simple  algebraic  stress  closures.  It  is  noteworthy  that,  except  for  the 
full  Reynolds  stress  closure  predictions  by  Reece  [1976],  all  numerical 
studies  of  the  problem  have  neglected  convection  and  diffusion  contribu¬ 
tions  to  the  balance  of  u.u . ,  thus  implying  the  condition  of  a  homogen- 
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eous  turbulent  flow. 


1.4  Outline  of  the  Thesis 


The  remainder  of  this  thesis  describes  the  study  performed  in 
greater  detail  throughout  the  next  six  chapters.  In  Chapter  2  the 
experimental  program  is  presented.  This  includes  descriptions  of  the 
apparatus  and  instrumentation,  the  experimental  methodology,  and  the 
major  error  sources  with  estimates  of  their  magnitudes.  In  Chapter  3 
the  numerical  model  framework  for  predicting  turbulent  flows  is  estab¬ 
lished.  In  Chapter  4  the  concept  of  partially-parabolic  flows  is  intro 
duced,  followed  by  a  detailed  account  of  the  numerical  procedure  as 
well  as  the  boundary  conditions  used. 

The  performance  of  the  numerical  procedure  and  the  turbulence 
models  were  separately  checked,  by  reference  to  calculations  of  several 
appropriately  documented  flows.  The  results  of  these  tests  are  pre¬ 
sented  in  Chapter  5.  Chapter  6  discusses  the  measured  and  calculated 
results.  Finally,  major  conclusions  and  specific  recommendations  are 
summarized  in  Chapter  7. 


CHAPTER  2.  THE  EXPERIMENT 


2.1  Introduction 

In  this  chapter  the  experimental  part  of  the  research  project  is 
described.  Section  2.2  outlines  the  flow  system,  the  instrumentation, 
and  the  flow  conditions  for  which  data  were  taken.  Section  2.3  details 
the  experimental  methodology  and  finally.  Section  2.4  discusses  major 
error  sources  and  estimates  of  their  magnitudes. 

2.2  Flow  System  and  Instrumentation 

The  experimental  system  was  composed  of:  a  water  rig,  of  which  the 
most  important  component  was  the  flow  test  section;  a  laser-Doppler 
velocimeter  and  its  associated  electronic  instrumentation;  a  motorized 
traversing  mechanism;  and  a  PDP-ll/34a  Digital  Equipment  Corporation 
minicomputer. 

The  basic  components  of  the  flow  test  section  and  the  coordinate 
systems  are  shown  schematically  in  Figure  2.1,  and  comprised  two  straight 
ducts  and  a  bend  of  square  cross-sections.  The  tangents  were  each  31 
hydraulic  diameters  long  and  were  respectively  attached  to  the  0°  (inlet) 
and  180°  (outlet)  planes  of  the  bend.  The  ratio  of  bend  mean  radius  of 
curvature  to  hydraulic  diameter  was  RC/DH  =  3.35.  The  bend  component 
was  constructed  by  machining  an  open,  curved  channel  of  mean  radius  of 
curvature  Rc  =  14.92  cm  (+0.02  cm)  into  one  of  the  faces  of  a  large, 
flat,  solid  piece  of  plexiglass  7.6  cm  thick.  A  plexiglass  plate  1.27  cm 


10 


rr 


thick  was  machined  to  fit  snugly  as  a  lid  over  the  open  channel  piece, 
so  yielding  an  enclosed  curved  duct  shape  of  cross-section  dimensions 
4.45  x  4.45  cm  (+0.02  x  0.02  cm  ).  An  0-ring  seal  placed  between  the 
open  channel  piece  and  the  lid  prevented  water  leaks.  The  modular  form 
of  this  construction  made  disassembly  easy  for  surface  cleaning  pur¬ 
poses.  The  upstream  and  downstream  tangents  were  each  138  cm  (+0.02 
cm)  long  and  were  also  constructed  from  plexiglass,  with  flat  walls 
1.27  cm  thick.  These  two  ducts  were  joined  by  flanges  to  the  bend, 
with  special  care  taken  to  avoid  possible  mismatches  between  the  com¬ 
ponent  cross-sections  which  otherwise  might  disturb  the  flow. 

A  flow-straightening  section  7.15  hydraulic  diameters  long  was 
placed  upstream  of  the  straight  duct  attached  to  the  bend  inlet  plane. 
The  purpose  of  this  device  was  to  uniformize  and  accelerate  the  devel¬ 
opment  of  the  cross-stream  plane  distribution  of  the  flow  approaching 
the  bend  in  the  upstream  tangent.  Various  arrangements  of  differently 
sized  stainless  steel  screens  were  tested  in  combination  with  one  or 
more  perforated  plexiglass  plates  3.175  cm  thick,  containing  85  holes  of 
3.175  mm  diameter  arranged  in  a  rectangular  array  spaced  4.495  mm  on  the 
centers  in  each  direction.  The  most  successful  arrangement  of  plates 
and  screens  in  the  uni formi zing  section  was  found  experimentally  and  is 
shown  in  Figure  2.2. 

The  test  section  was  part  of  a  closed  loop  system  through  which 
water  at  20°C  was  made  to  flow  by  gravity  from  a  constant  head  tank. 

From  this  tank  the  flow  passed  through  the  test  section,  and  then  into 
a  large  sump  tank  from  where  it  was  pumped  back  to  the  constant  head 
tank  by  a  3/4  HP  Burkes  centrifugal  pump.  The  constant  head  was  ensured 
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by  a  large  diameter  PVC  overflow  pipe  extending  through  the  bottom  of 
the  head  tank.  Flow  to  the  head  tank  was  controlled  by  a  gate  valve 
and  measured  using  a  2  in.  Barco  venturi  meter  (P.N.  2-393)  connected 
to  a  50  in.  (1.27  m)  differential  mercury  manometer.  In  order  to  rule 
out  the  possibility  of  propagating  perturbations  induced  by  flow  com¬ 
ponents,  the  use  of  valves,  sharp  bends  and  metering  devices  was  avoided 
altogether  along  the  test  section  flow  loop.  Flow  to  and  from  the  test 
section  tangents  was  channelled  through  2  in.  i.d.  tygon  tube  pieces; 
flexible  enough  to  be  bent  without  kinks  over  a  large  radius  of  curva¬ 
ture,  yet  stiff  enough  to  avoid  wall  collapse  due  to  flow-induced  pres¬ 
sure  drop.  Baffles  located  in  the  constant  head  tank  served  to  dampen 
the  swirling  motion  of  the  flow  leaving  the  tank.  Residual  swirl  from 
the  head  tank  and  secondary  motions  induced  by  mild  curvature  in  the 
tygon  tube  upstream  of  the  test  section  were  eliminated  by  the  flow- 
uniformizing  section  placed  between  the  tygon  tube  and  the  upstream  tan¬ 
gent.  All  experiments  were  conducted  for  the  flow  rate  condition  imposed 
by  the  constant  head  tank.  This  corresponded  to  a  Reynolds  nun*^**  of 
Re  *  56,700  and  a  Dean  number  of  De  *  21,900  in  the  flow  test  section. 

Measurements  of  the  mean  flow  and  turbulence  characteristics  were 
made  using  the  laser-Doppler  velocimeter  technique  in  backscatter  mode. 
The  velocimeter  employed  is  shown  schematically  in  Figure  2.3  in  rela¬ 
tion  to  the  flow  test  section.  It  comprised  a  2  watt  Lexel  Argon-Ion 
water-cooled  laser,  a  mirror  stage  for  reflecting  the  laser  beam  180° 
into  the  velocimeter  optics,  the  optics,  and  a  4  in.  (10.2  cm)  diameter 
mirror  for  reflecting  the  converging  velocimeter  beams  from  the  horizon¬ 
tal  to  the  vertical  direction.  This  mirror  also  served  the  function  of 
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reflecting  the  back-scattered  Doppler-shifted  radiation  into  the  veloci- 
meter  collecting  optics.  The  optics  were  of  the  DISA  55X  Modular  Series 
and  consisted  of:  two  separately  adjustable  quarter-wave  retardation 
plates;  a  50:50  neutral  beam  splitter;  a  beam  color  splitter;  a  back- 
scatter  unit  (containing  a  mirror  inclined  45°  with  respect  to  the 
velocimeter  optical  axis,  and  serving  to  support  at  right  angles  to  the 
optical  axis  the  photomultiplier  optics  consisting  of:  a  color  separa¬ 
tor;  two  interference  filters  and  two  RCA-4526  photomultiplier  tubes); 
a  pinhole  section;  a  beam  translator  and  a  310  mm  achromatic  focusing 
lens. 

The  laser  and  velocimeter  optics  were  mounted  to  the  top  of  a  thick 
aluminum  table,  which  was  itself  firmly  bolted  to  an  x,  y,  z  traversing 
mechanism.  The  traversing  mechanism  could  displace  the  table  top  +  7.5 
cm  in  5  urn  increments  along  any  of  the  coordinate  axes  by  means  of  three 
linearly  encoded  stepping  motors  monitored  by  the  Digital  Equipment 
Corporation  PDP  11 /34a  minicomputer.  The  minicomputer  functioned  as 
the  central  data  acquisition  and  reduction  controller.  In  addition  to 
directing  the  spatial  sequence  of  an  experimental  run,  the  computer  was 
programmed  to  conduct  the  acquisition,  statistical  processing,  plotting 
and  storage  of  Doppler  data  validated  and  measured  by  a  DISA  55L96 
Doppler  signal  processor  or  "counter".  The  PDP  11/ 34a  has  a  256  K  16 
bit  random  access  memory  and  is  equipped  with  dual  hard  RL01  magnetic 
disc  drives  (5  Mbytes  each).  The  computer  interacts  by  means  of  an  RT-11 
software  package  with  various  input-output  devices,  including  a  Tektroni 
4025  graphics  terminal,  a  Decwriter  II  hardcopy  terminal  and  a  Tektronix 
4662  digital  pen  plotter. 
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2.3  Experimental  Methodology 

Prior  to  an  experimental  run,  water  was  allowed  to  flow  through  the 
rig  until  it  was  purged  of  air  bubbles  and  had  attained  a  steady  thermal 
state  corresponding  to  20°C  {+  1°C).  Mass  flow  through  the  test  section 
was  controlled  by  setting  the  constant  head  overflow  condition  to  a  mere 
trickle  and  continuously  monitoring  the  pressure  drop  through  the  ven¬ 
turi  meter  connected  to  the  head  tank  feed  line. 

At  any  given  streamwise  measurement  station,  the  velocimeter  sup¬ 
port  table  was  manually  positioned  such  that  the  velocimeter  horizontal 
optical  axis  was  oriented  perpendicular  to  the  test  section  side;  see 
Figure  2.3  for  an  example  corresponding  to  a  streamwise  location  of 
e  =  180°.  Fine  adjustments  to  the  90°  beam  deflector  mirror  ensured 
that  the  velocimeter  vertical  optical  axis  was  perpendicular  to  the  test 
section  top  surface.  In  combination,  these  adjustments  produced  a  fringe 
pattern  which  was  parallel  (to  within  +0.3°  at  any  y,  z  cross-section 
location)  to  the  cross-stream  plane  of  the  test  section.  The  velocimeter 
optical  probe  volume  was  formed  by  the  intersection  of  two  514.5  nm 
(green)  light  beans  with  a  half-angle  in  air  of  4.90°,  for  which  the 
volume  characteristics  were:  a  diameter  of  0.09  mm,  a  length  of  1.1  mm, 
and  a  fringe  spacing  of  3.02  ym  with  about  28  fringes  contained  in  the 
probe.  In  reality,  spatial  filtering  and  threshhold  settings  on  the 
counter  reduced  the  dimensions  of  the  optical  probe.  The  probe  volume 
was  positioned  at  the  top  outer-wall  corner,  inside  the  test  section, 
by  fine  control  of  the  motorized  traversing  table.  Positioning  in  the 
x  and  y  coordinate  directions  was  accurate  to  within  +0.05  mm  while 
positioning  in  the  z  direction  was  accurate  to  better  than  +0.5  mm. 


With  the  reference  corner  position  at  a  streamwise  location  estab¬ 
lished,  the  computer  software  was  activated  which  controls  signal  acqui¬ 
sition  and  data  processing  on  a  sequentially  scanned  measurement-grid. 
After  flow  symmetry  had  been  established  at  various  strearrwise  stations, 
the  bulk  of  the  measurements  were  restricted  to  a  symmetrical  half  of 
the  test  section,  on  a  grid  consisting  of  4  to  5  profiles  at  different 
z  locations,  each  containing  29  to  31  equally  spaced  points  in  the  y 
direction.  At  each  point  on  the  measurement-grid  the  mean  flow  and 
turbulence  characteristics  were  statistically  determined  from  popula¬ 
tions  of  5  to  10  samples  consisting  of  1,000  measurements  each.  Each 
measurement  was  required  to  satisfy  the  counter  5/8  validation  compari¬ 
son  to  within  a  preset  tolerance  of  3%.  At  every  validation  of  a  Doppler 
burst  a  "data  ready"  signal  was  issued  by  the  counter  to  a  logic  conver¬ 
sion  circuit.  This  circuit  then  sent  a  triggering  pulse  to  the  computer 
parallel  line  interface  module  which  was  checked  for  data  availability 
by  a  software  loop  approximately  every  20  ys.  No  interrupt  routines 
were  used  for  obtaining  data  due  to  the  higher  sample  rates  made  possi¬ 
ble  by  the  handshake  technique.  Data  rates  of  about  1  kHz,  with  approxi¬ 
mately  60%  validation,  were  obtained  after  the  flow  was  seeded  with  corn¬ 
starch  particles  ranging  in  size  between  1  and  10  ym.  Optical  alignment 
and  automatic  grid  scanning  were  performed  at  the  following  streamwise 
locations:  XH  *  -5,  -1,  1 ,  5,  10,  20  in  the  straight  ducts,  and  6=3°, 
45°,  90°,  130°,  177°  in  the  bend. 

Although  capable  of  two-component  measurements;  the  availability  of 
only  one  counter  during  the  early  stages  of  this  work  restricted  the  use 
of  the  velocimeter  to  single  channel  mode  for  most  of  the  experiments. 


Values  for  the  streamwise  component  of  mean  velocity  and  normal  stress, 
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UQ  and  Ug ,  were  derived  directly  from  measurements  obtained  with  the 

velocimeter  fringes  aligned  perpendicularly  to  the  streamwise  coordinate 
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direction.  Values  for  the  transverse  components,  Ur  and  ur,  and  for  the 

Reynolds  stress,  u_u.  were  derived  as  described  in  Melling  and  Whitelaw 

[1976]  by  combining  measurements  obtained  with  fringes  oriented  at  angles 

of  +  45°  and  -  45°  with  respect  to  the  streamwise  direction.  As  a  check, 

-y 

measurements  of  UQ  and  uQ  obtained  in  this  manner  were  found  to  agree  to 
within  experimental  error  with  the  direct  observations.  Measurements 
involving  the  velocity  component  in  the  spanwise  (z)  coordinate  direc¬ 
tion  could  not  be  made  accurately  due  to  optical  inaccessibility  of  the 
flow. 

While  the  bulk  of  the  measurements  were  made  in  the  manner  described 
above,  with  the  velocimeter  in  single  channel  mode,  towards  the  end  of 
the  experiment  the  availability  of  a  second  counter  allowed  us  to  per¬ 
form  a  limited  number  of  measurements  using  the  two  velocimeter  channels 
simultaneously.  For  this  case  the  interference  fringe  patterns  were 
respectively  aligned  parallel  and  perpendicular  to  the  streamwise  velo¬ 
city  component.  To  resolve  flow  directional  ambiguity,  and  also  to 
optimize  the  filter  range  of  the  counters,  a  net  frequency  shift  of  700 
kHz  was  imposed  on  both  channels  using  a  DISA  55N10  Brag  cell  combined 
with  electronic  downmixing.  In  this  way  additional  measurements  at 
XH  =  -5,  -1,  1  and  5  in  the  respective  tangents  were  obtained.  In  addi¬ 
tion,  careful  checks  were  conducted  of  the  data  previously  obtained  at 
XH  =  -5  and  5  and  at  6  =  3°  and  177°.  In  general,  the  checks  showed 
very  good  agreement  between  the  two  measurement  methods,  for  both  of  the 


velocity  components  and  their  respective  normal  stresses.  However,  at 
=  -5  the  two  methods  differed  markedly  with  respect  to  the  measure¬ 
ment  of  the  very  weak  transverse  velocity  component  arising  at  this 
location.  The  first  method  consistently  showed  larger  scatter  in  the 
measurements  and  the  uncertainty  was  ultimately  traced  to  relatively 
small  but  significant  inaccuracies  associated  with  determining  the 
+45°  orientations  required  by  this  method.  At  XH  =  -5,  -1,  1  and  5, 
and  6  =  3°,  it  is  the  more  accurate  frequency-shifted  transverse 
velocity  data  which  is  reported  here. 

2.4  Error  Estimates 

Error  sources  affecting  the  accuracy  (systematic  error)  and  preci¬ 
sion  (random  error)  of  laser-Doppler  measurements  have  been  discussed 
by,  for  example.  Durst,  Mel  ling  and  Whitelaw  [1976],  Drain  [1980]  and 
Buchhave  [1979].  In  this  study  the  most  serious  systematic  errors  were 
attributed  to  velocity  gradient  broadening  and  velocity  bias  respective¬ 
ly.  Velocity  gradient  broadening  has  been  analyzed  by  Melling  [1975] 
who  proposed  a  simple  method  for  estimating  its  magnitude.  Various 
weighting  methods  have  been  proposed  by,  among  others,  McLaughlin  and 
Tiederman  [1973],  George  [1975]  and  Buchave  [1975]  to  correct  for  the 
velocity  bias  effect,  but  none  of  these  is  entirely  satisfactory;  they 
all  involve  assumptions  regarding  the  statistical  distribution  of  par¬ 
ticles  in  the  flow  and,  in  practice,  the  corrections  can  be  influenced 
by  the  additional  problem  of  "Incomplete-signal  bias".  For  the  condi¬ 
tions  of  this  study,  gradient  broadening  and  velocity  bias  were  estimated 
to  be  significant  only  in  the  near  wall  regions  of  the  flow,  where 
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velocity  gradients  and  turbulence  intensities  were  highest  and  the  data 
rate  lowest.  Fortunately,  the  errors  are  of  similar  magnitude  and  of 
opposite  sign,  tending  to  cancel  out  their  respective  effects  on  the 
measurements.  For  this  reason,  and  because  the  errors  were  small  anyway, 
corrections  were  not  applied  to  the  measurements.  Table  2.1  presents 
estimates  of  the  maximum  combined  inaccuracies  of  these  two  error 
sources  on  the  quantities  measured. 

Values  of  the  transverse  mean  velocity  component,  Ur>  and  of  the 
Reynolds  stress,  u^u”,  were  prone  to  a  third  systematic  error.  At  the 
last  four  locations  in  the  bend,  and  at  XH  *  10  and  20,  these  two 
quantities  were  derived  from  measurements  taken  at  +  45°  and  -  45°  to 
the  streamwise  flow  direction.  An  error  in  setting  this  90°  angle 
could  seriously  affect  the  accuracy  of  these  measurements.  The  problem 
has  been  considered  by  Humphrey  [1977]  for  a  flow  of  similar  character¬ 
istics  to  the  present  one,  and  he  shows  that  an  angular  uncertainty  of 
0.4°  can  lead  to  an  error  of  5%  in  U  and  +  3%  in  uQu  .  In  this  work 
special  care  was  taken  to  ensure  a  maximum  angular  uncertainty  of  less 
than  0.3°  in  setting  the  +  45°  bean  orientations. 

The  two  main  sources  of  random  error  affecting  the  precision  of 
the  measurements  were  attributed  to  statistical  sampling  uncertainty 
(due  to  the  finite  size  of  sample  populations)  and  uncertainty  in  the 
determination  of  the  reference  or  normalizing  velocity,  Ug.  Estimates 
of  the  first  uncertainty  were  derived  from  the  measurements  themselves 
and,  for  all  quantities,  were  found  to  be  less  than  +  1%  r.m.s.  error. 

The  error  in  Ug  was  larger  {+  2%)  and  arose  principally  from  uncertain¬ 
ties  in  the  construction  of  the  venturi  meter.  As  a  further  check,  the 


bulk  mass  flow  at  each  longitudinal  station  was  estimated  by  integrating 
the  measurements.  This  yielded  a  value  of  Ug  which  was  within  +_ 6 % 
of  the  venturi  meter  measurement.  Estimates  of  the  maximum  combined 
effects  of  these  two  errors  on  the  measurements  are  provided  in  Table 


Measurements  of  the  pressure  coefficient,  Cp,  in  the  straight  and 
curved  duct  test  sections,  using  side  wall  pressure  tappings  connected 
to  an  inclined  manometer  bank,  were  also  prone  to  a  random  r.m.s.  error 
ranging  from  +  10%  at  low  absolute  values  of  Cp  to  +  5%  at  the  higher 
values.  The  uncertainty  in  Cp  was  due  mainly  to  the  reading  error 
associated  with  the  manometer  bank. 


CHAPTER  3.  MEAN  FLOW  EQUATIONS  AND  TURBULENCE  MODELS 


3.1  Introduction 

In  this  chapter  the  numerical  modeling  framework  for  predicting 
turbulent  flows  is  detailed.  Section  3.2  presents  the  mean  flow  equa¬ 
tions  and  discusses  the  problem  of  closure  common  to  all  nonlinear  sto¬ 
chastic  systems.  Section  3.3  describes  the  models  tested  in  the  present 
study.  These  include  the  k-e  model  of  turbulence  in  Section  3.3.1  and 
the  algebraic  stress  model  (ASM)  in  Section  3.3.2.  The  model  constants 
for  both  the  k-e  and  the  ASM  closures  are  also  given  in  Section  3.3.2. 
Finally,  the  ASM  relations  appropriate  for  the  present  flow  configuration 
is  summarized  in  Section  3.3.3. 

3.2  Mean  Flow  Equations  and  the  Problem  of  Closure 

In  the  present  study  numerical  calculations  of  turbulent  flows  were 
based  upon  the  time-averaged  Navier-Stokes  equation  and  continuity  equa¬ 
tion  as  first  proposed  by  Osborne  Reynolds.  For  a  statistically  station¬ 
ary  flow  of  a  fluid  with  uniform  density  p,  Reynolds  decomposition  of  the 
field  variables  followed  by  time  averaging  of  the  resulting  equations 
yields 

Continuity; 


(3.2.1) 


Momentum: 
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where  lower  and  upper  case  u's  stand  for  fluctuating  and  time-averaged 
velocity  components  respectively,  P  represents  mean  pressure,  overbars 
imply  time  averaging  of  the  correlations  in  question  and  y  denotes  the 
viscosity  of  the  fluid. 

Equations -(3.2.1 )  to  (3.2.4)  can  be  applied  to  both  cylindrical 
(r,  0,  z)  and  rectangular  (y,  x,  z)  coordinate  systems.  To  get  the 
appropriate  equations  in  rectangular  coordinates,  set  j  =  0  and  r  =  1 
and  make  the  following  substitutions: 


3? “"‘ax  ’  Jr  ~*“Jy  *  Ue"fc  Ux  and  Ur*~  Uy 

wherever  they  appear.  The  governing  equations  in  cylindrical  coordinates 

can  be  deduced  simply  by  setting  j  =  1  in  equations  (3.2.1)  to  (3.2.4). 

Unfortunately,  the  above  set  of  equations  can  not  be  solved  directly 

for  the  mean  velocities  and  pressure  due  to  the  appearances  of  the  six 

correlations:  u  u  ,  uTuT,  u  u  ,  ZTu7,  u  u  and  uQu  .  The  quantities 
rr  t?  u  z  z  ro  rz  oz 


-pu  u  ,  -Pu-u.,  -pu  u  ,  -PU  uQ,  -pu  u  and  -pu_u„  are  known  as  turbulent 
or  Reynolds  stresses.  Exact  transport  equations  for  these  correlations 
can  be  derived,  however  they  contain  correlations  of  even  higher  order 
due  to  the  nonlinear  nature  of  the  equations.  "Closure"  therefore  can 
not  be  achieved  by  resorting  to  solving  transport  equations  of  higher 
and  higher  order.  Instead,  model  approximations  must  be  introduced  at 
a  certain  order  in  terms  of  lower  order  correlations  and  mean  quantities. 
This  constitutes  the  task  of  turbulence  modeling  which  is  the  subject  of 
next  section. 
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3.3  Modeling  of  Turbulence 

As  mentioned  earlier,  the  role  of  turbulence  modeling  is  to  provide 
a  path  for  the  determination  of  the  Reynolds  stresses  appearing  in  equa¬ 
tions  (3.2.2)  to  (3.2.4).  Depending  upon  the  level  of  closure  and  the 
generality/complexity  desired,  a  wide  variety  of  models  have  been  pro¬ 
posed  and  tested.  For  a  comprehensive  review,  see  Launder  and  Spalding 
[1972],  Reynolds  [1976]  and  Rodi  [1978].  Two  models  of  turbulence  were 
selected  for  the  present  study:  a  two-equation  (k-e)  model  and  an 
algebraic  Reynolds  stress  model.  Their  essential  features  are  outlined 
in  the  next  two  subsections. 
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3.3.1  k-e  Model 

Following  Boussinesq  [1877],  the  turbulent  stresses  are  related  to 
the  mean  rate-of-strain  tensor  via  the  definition  of  an  isotropic  tur¬ 
bulent  viscosity,  y^.  In  cylindrical/rectangular  coordinates,  these 
relations  take  the  forms  given  below. 


-fWr-  M2!?)  "  ffk  (3.3.1. 1) 
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-popix  -  W  +4t)  (3-3-1-6) 

Equations  (3. 3. 1.1)  to  (3. 3. 1.6)  are  completely  analogous  to  the 
constitutive  relations  for  an  incompressible  flow  of  a  Newtonian  fluid, 
except  for  the  additional  term  2pk/3  which  appears  in  each  of  the  normal 
stress  relations.  This  addition  is  necessary  to  ensure  that  the  defini¬ 
tion  of  turbulence  energy 

k  -  i(u  u  +  u„u„  +  u  u  ) 

2  r  r  0  0  z  z 

is  not  violated. 

Substituing  equations  (3. 3. 1.1)  to  (3. 3. 1.6)  into  equations  (3.2.2) 
to  (3.2.4)  and  rearranging  yields: 


H  (r^O  iwi  ♦  &  (f aw  -S*#  -  - 


(3. 3.1. 7) 


_l£' 

rig 

■2^)} 

j[2®‘ft^-£(pt4>«*(^*-JA)] 

(3. 3. 1.8) 

!-^(frum)«^k(fuoI)*^(fifcu,)  -  -|f' 

1±£im  &>*£(* 

(3.3.1. 9) 

9'  -  P  +  |fk 

Ht  -  Mt  *  b  . 

With  the  Boussinesq  eddy  viscosity  approximation,  the  task  of 


modeling  is  shifted  to  the  determination  of  y^.  Using  the  k-e  model, 
ut  is  determined  from  two  turbulence  quantities:  the  turbulence  energy 
k  and  its  rate  of  dissipation  e,  via  the  relation 


where  Cy,  to  first  approximation,  is  a  constant  of  proportionality. 
Following  Jones  and  Launder  [1972],  the  modelled  k  and  e  transport 
equations  are: 
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Equations  (3.2.1)  and  (3. 3. 1.7)  to  (3.3.1.13),  together  with  appro¬ 
priate  boundary  conditions  form  a  closed  system.  Results  of  k-e  model 
calculations  will  be  presented  in  Chapter  6. 

3.3.2  Algebraic  Stress  Model 

Algebraic  stress  models  (ASM's)  are  special  cases  of  full  Reynolds 
stress  models  proposed  by,  among  others,  Hanjalic  and  Launder  [1972] 
and  Launder  et  al.  [1975]. 

The  starting  point  of  ASM  is  the  transport  equation  for  u.u..*  For 

*  J 

a  fluid  of  uniform  density  and  viscosity  and  unaffected  by  external  force 
fields,  this  equation  can  be  written  as: 


*For  ease  of  presentation  and  general  discussion,  Cartesian  tensor  nota¬ 
tion  is  adopted  in  this  section.  Final  algebraic  relations  appropriate 
to  boch  cylindrical  and  rectangular  coordinates  are  summarized  in 
Section  3.3.3. 
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Equation  (3.3.2.T)  is  obtained  by  multiplying  the  instantaneous 
x.  momentum  equation  by  u.  and  adding  it  to  the  instantaneous  x.  momen- 
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turn  equation  multiplied  by  u.  then  time-averaging  the  resultant  equa¬ 
tion.  The  physical  process  represented  by  each  term  is  also  indicated, 
for  a  detailed  explanation,  see  Tennekes  and  Lumley  [1972]. 

With  closure  at  the  second-order  (u,u .)  level,  the  production  term 

'  J 

contains  only  turbulent  stresses  and  mean  velocity  gradients  and,  there¬ 
fore,  requires  no  modelling.  The  dissipation  and  pressure-strain  correla¬ 


tion  terms  were  modelled  as  in  Gibson  and  Launder  [1978]  with 


Ci3  =  t  5ij  E 


(3. 3. 2. 2) 


*ij  *  ^ij  ,1  +  4ij,2  +  ^i j ,1  +  ^i j ,2 


(3. 3. 2. 3) 


Equation  (3. 3. 2. 2),  first  proposed  by  Rotta  [1951],  is  based  upon 


the  supposition  that  the  dissipative  motions  are  isotropic  for  flows  at 
relatively  high  Reynolds  numbers.  The  pressure-strain  term  plays  an 
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influential  role  in  the  u.u.  equations.  Since  this  term  makes  no  net 

contribution  to  the  turbulence  energy  but  merely  redistributes  it  among 

the  normal  stresses,  it  is  also  referred  to  as  the  "redistribution" 

term.  The  theoretical  background  for  modelling  <p .  ■  in  terms  of  equa- 

i  J 

tion  (3. 3. 2. 3)  was  discussed  in  Hanjalic  and  Launder  [1972]  and  Launder 

et  al.  [1975].  It  can  be  shown  that  there  are  three  types  of  physical 

processes  which  dictate  the  distribution  of  <)>..:  1)  mutual  interactions 

^  J 

between  turbulence  components  (4>. .  ,);  2)  mean  rate-of-strain  interacting 

•  J  9  * 

with  the  turbulence  (*. .  9);  and  3)  corrections  to  <f>..  ,  and  <j>..  9 
resulting  from  the  presence  of  walls  ($'..  ■,  and  <j>'..  ,).  The  modelling 
of  the  four  terms  in  equation  (3. 3.2. 3)  is  outlined  next. 

For  *..  , ,  Rotta's  [1951]  proposal  was  adopted  whereby 

1  J  »  I 

*«.)  ■  k>  •  '3-3-2-4> 
The  primary  functions  of  this  term  are  to  equalize  normal  stresses  and 
to  diminish  to  zero  shear  stresses.  For  this  reason,  it  is  usually 
referred  to  as  the  "return- to-isotropy"  term. 

For  <j>..  Launder  et  al.  [1975]  have  suggested 

1  J  9  ^ 

*1j.2  *  -C2(Pij  -  I  P>  (3-3-2'5’ 

as  the  simpler  alternatives  for  a  more  complete  model  of  which  the  first 
term,  apart  from  a  proportionality  constant,  is  the  same  as  that  given 
in  equation  (3. 3. 2. 5). 

Experiments  show  that  the  presence  of  a  wall  affects  the  turbulence 
by  damping  the  fluctuating  velocity  component  normal  to  the  wall  and  by 
enhancing  those  components  along  the  other  two  directions.  In  the 
present  study  these  wall  effects  were  modelled  through  ,  and  <J>/.  9  as 

I  J  9  *  I  J  >  ^ 
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4ij,l  =  Ci(l>  (Vm  Vm  {ij  ‘  f  Vi  Vj  *  !  Vj  Vi)f(F77) 


’ij,2  *  C2(*km  Vm  Sij  "  2  *ki  nk"j  '  2  *kj  Vi1  f(n*rj* 


(3. 3. 2. 6) 


(3. 3. 2. 7) 


where  r.  is  the  position  vector  to  the  point  in  question,  £  is  a  charac¬ 
teristic  turbulent  length  scale,  n^  is  the  unit  normal  to  the  wall  and  f 
is  an  empirical  function  of  distance  whose  primary  role  is  to  diminish 
the  respective  influences  of  q1..  ,  and  q'..  7  with  increased  distance 
from  the  wall.  Equation  (3. 3. 2. 6)  was  first  proposed  by  Shir  [1973]  in  a 
numerical  study  of  atmospheric  turbulent  flows  in  the  idealized  planetary 
boundary  layer.  Gibson  and  Launder  [1978]  later  extended  Shir's  idea  to 
account  for  the  wall  effect  expressed  by  q'. .  9  as  given  by  equation 

(3. 3. 2. 7).  For  a  two-dimensional  straight  channel,  f  was  assumed  to 
vary  as  recommended  by  Reece  [1977]  according  to  the  expression 

f(7)  =  ^r[7  +  ^  (3-3-2-8) 


where  y  is  the  distance  measured  from  one  of  the  two  walls  and  D  is  the 


channel  width.  The  constant  Cw  was  chosen  such  that  f  -*■  1  as  y  -*■  0.  It 

can  be  shown  that  in  a  region  where  the  logarithmic  law  of  the  wall  pre- 

3/4 

vails  and  the  turbulence  is  in  local  equilibrium  that  Cw  *  ic/Cy  .  For 
flows  inside  a  square  duct  with  two  sets  of  opposing  walls,  it  was 
assumed,  following  Reece  [1977],  that  the  influences  of  both  sets  of 
walls  on  4>..j  are  independent  of  each  other  and  can  be  added  algebraically. 

The  model  assumptions  introduced  so  far  in  equations  (3. 3.2.2)  to 

(3. 3. 2. 8)  for  and  are  algebraic  since  they  contain  no  spatial 


gradient  of  the  Reynolds  stresses.  If  the  remaining  convection  and 


diffusion  terms  in  the  u.u  .-transport  equation  can  also  be  modelled  by 

I  J 

expressions  without  gradients  of  u.u . ,  equation  (3. 3. 2.1)  is  reduced 
to  the  following  algebraic  form 


Vj  =  uiuJ(Vq 


»  k  i  c ) . 


(3. 3. 2. 9) 


Models  of  this  kind  have  been  proposed  by  Launder  [1972]  and  Rodi  [1976] 
and  it  was  the  more  general  model  of  Rodi  which  was  adopted  in  the 


present  study. 

Following  Rodi  [1976],  the  convection  and  the  diffusion  terms  were 
modelled  collectively  as: 


D  u.u.  _ _  u.u.  n.  u.u. 

-tH  -  ■  -Y  $  -  w»>]  *  -¥  (p  -  <3-3-2->°> 

where  &{  )  stands  for  "diffusion  of  the  quantity  in  parenthesis."  The 
second  part  of  equation  (3.3.2.10)  results  from  mere  definition.  Sub¬ 
stituting  all  the  modelling  assumptions  into  equation  (3. 3.2.1)  and 
rearranging  yields 


.Vi 


"  3  5ij  =  e  {(Pij  "  3  6ij  P)  +  1-C9  (*’ij,l  +  *'ij,2)} 


1  -  C, 


(3.3.2.11) 


p-v7£ 


and  d>'  .  , 


and  9  given  by  equations  (3. 3. 2. 6)  to  (3. 3. 2. 8). 


It  is  also  necessary  to  provide  a  means  for  the  determination  of 
k  and  e  and  this  was  achieved  by  solving  modelled  transport  equations 


3 


for  both  variables.  Following  Launder  et  al.  [1975],  the  k  and  e 
equations  adopted,  in  cylindrical/rectangular  coordinates,  are: 


i £ (frUrkj  +  &  (f  u.  k ) «  £  (f  ife Id  -  f  p  -  f  £. 

♦  ♦  Sk  (3.3.2.12) 

-Cup-I-1 


Equations  (3.3.2.12)  and  (3.3.2.13)  are  identical  to  those  used  for  the 
k-e  model,  except  for  diffusion  terms  which  are  now  approximated  without 


introducing  the  notion  of  an  isotropic  turbulent  viscosity,  u^. 

The  determination  of  the  model  constants  appeared  in  Sections 
3.3.1  and  3.3.2  have  been  made  by:  1)  available  experimental  data  for 
simple  flows;  and  2)  computer  optimization.  For  a  comprehensive  review, 
see  the  lecture  notes  of  Launder  [1980].  The  set  of  model  constants 
adopted  for  the  present  study  is  listed  below: 


<  - 
C  : 


0.4187 
•  0.09 


a,  =1.0 


a  = 


<2/(C 


9  -  c  ,)/C 

e2  el  y 


1/2 


'el 


=  1.44 


Ce2=1.92 


'e3 


kl 

C,  = 


■  0.36  (Ce2  -  C£l) 

•  0.22 
1.8 


C0  =  0.6 


r  1 

h 

C  1 


0.5 

0.3 


3.3.3  Summary  of  the  ASM  Relations  for  Flows  in  a  Duct  of  Square 
Cross-Section 

The  ASM  outlined  in  Section  3.3.2  with  Cartesian  tensor  notation 
is  summarized  here  for  both  the  straight  and  the  curved  flow  configura 
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tions  of  interest  to  this  study. 

The  expressions  for  e. P..  ,  and  ,,  when  expanded  in 

cyl indrical/rectangular  coordinates,  take  the  forms  given  below: 


c,'  (7)  juT^F  -2uT^  <k  ] 
c!  (7)  |  uto;  (k +  iMx  f  ] 
c!  (7)  I  Q  -2  uTuzf] 
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Qri  = 

C.'  l£) 

[-1.5  u,u,  (.F 

o] 

<  & 

[-L5  k»uj  Fj 

<*" 

-  CtC i  1 

(P,*-fP)F 

“  2  (  Pnr 

-fw«J 

-CjCj  ( 

lPrr-|p)Q 

♦CRr* 

|p)f} 

<2  - 

-CtCt‘  { 

(Prr-fp)<i 

-2(Pn 

■fWF} 

- 

1.5  CiCt  Prt  C, 

(firz.J  ~ 

1.5  CiClPn  IF+G) 

L5  CiC,'P„F 

where 


k3/2/i 


U-sWl 


k»'*/e  r  i 

c,  l  y 


U  *-i_l 
1  y  o«-y  J 


y  «  (r-rj 


Straight-forward  but  tedious  manipulation  of  the  above  expressions 


result  in  the  following  algebraic  relations  for  the  six  unknowns  err, 
eee-  ez2  ers ’  erz  and  eez: 


ft.  1-841^ 

^  ♦  2  U*2A2)  ♦  2A4  -A3 

r  - 1.5  M  J  +  ft»l  -Bfe(  ^  *J  ] 

£  ♦  L5A3]  4  M-66  ♦  B5l^-J  «*}] 

-Jt)-a3-A4-2A2  |£-2»l 

p  *I.3A4]+  ft.[B3l^*J^)] 
|^-3A3-iL]4ft.[K(^-j^)483^] 


2(I-*2A1) 


♦  2A3-A4 


« w  [  -B9  %#  ]♦  e«  f-B9(  %)]  *  €a  [  0  ] 

e„  [  B9  $•'  -  2.29  A4 ]♦  €„ [-69  J 

+  e"[-B7ls?  ]+  M'B7l^|£)-Z25(A3*M)-^] 
err[o]*CM[-8e^]^l[-B8^] 

♦e*[-68  &  ]  ♦  e,t-B8  ^  ]♦  Mae  ^r-“5 «- S  ] 

3  00  l  rtfi  az  ' 


Cl  cl 
I- Cl 


I- Cl  *  kJ  h 


tv*' 


t-c*  V 


61 

B2 

B3 

B4 

B5 

B& 

B7 

B8 


2Xk 
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2(U2Ai)  ♦  A2 
|  +  2(Ai  +  A2) 
1+2A1-A2 
2C1+2A2)**  A1 
I  ♦  2A2  -  A1 
2*Al«tA2 

1.5  (  1+  15A1  +  L5A2) 
1.5  (1-1.5  Al) 


*/• 


15  0*1.5  A2)  . 
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CHAPTER  4.  NUMERICAL  PROCEDURE 


4.1  Introduction 

In  view  of  the  complexity  of  the  governing  equations,  solutions  to 
the  problem  are  sought  via  numerical  methods.  In  the  present  study,  the 
structure  and  strategy  of  the  TEACH  calculation  procedure  developed  at 
Imperial  College  [1976]  are  adopted,  with  major  modifications  to  take 
advantages  of  the  partially-parabolic  nature  of  the  flow.  Section  4.2 
introduces  the  concept  of  partially-parabolic  flow  and  discusses  the 
consequence  of  this  classification  in  terms  of  the  describing  equations. 
Section  4.3  details  the  discretization  practices.  Section  4.4  summar¬ 
izes  the  boundary  conditions  for  the  calculations.  Finally,  Section 
4.5  outlines  the  solution  algorithm. 

4.2  Partially-Parabolic  Flows 

There  are  primarily  three  mechanisms  which  transmit  local  distur¬ 
bances  to  other  points  of  a  flow:  convection,  molecular  diffusion  and 
pressure  waves.  Strictly  speaking,  all  steady  subsonic  flows  are  ellip¬ 
tic  in  the  sense  that  perturbations  at  a  point  can  influence  the  flow 
state  at  any  other  point.  A  solution  to  this  type  of  problem  requires 
the  specification  of  boundary  conditions  on  a  closed  boundary  of  the 
flew  domain.  From  a  computational  point  of  view,  however,  two  addi¬ 
tional  classifications  are  possible:  parabolic  flows  and  partially- 
parabolic  flows. 


In  terms  of  the  three  mechanisms  mentioned  above,  a  flow  situation 
is  classified  as  being  partially-parabolic  if  along  a  certain  direction 
the  only  means  for  transmitting  local  disturbances  upstream  is  through 
the  pressure  field;  there  must  be  no  flow  reversal  and  diffusion  should 
be  negligible  along  the  "partially-parabolized"  direction. 

The  flow  configuration  of  interest  to  the  present  study  falls  into 
the  category  of  partially-parabolic  flows  as  long  as  the  bend  is  not  so 
tight  as  to  render  flow  reversal  along  the  main  flow  direction.  The 
problem  of  a  deflected  jet  or  a  jet  in  a  crosswind  flow,  as  studied  by 
Bergeles  [1976],  can  also  be  classified  as  pertaining  to  the  partially- 
parabolic  category. 

The  immediate  consequence  of  this  "partial-parabolization"  is  that 
all  the  underlined  terms  in  the  equations  presented  in  Chapter  3  are 
assumed  negligible  and  so  were  neglected  in  the  calculations.  The 
validity  of  this  simplification  was  checked  by  Humphrey  et  al .  [1981] 
for  turbulent  flow  in  a  square  duct  of  strong  curvature.  It  was  found 
from  their  calculation  that  the  longitudinal  diffusion  was  no  larger  than 
2%  of  the  longitudinal  convection.  It  would  appear,  therefore,  that  a 
partially-parabolic  procedure  accounting  for  strong  pressure  variations 
could  provide  more  precise  results  through  increased  grid  refinement. 

Since  all  transport  equations  given  in  Chapter  3  share  the  same 
general  form 

i  &  (frUr*) + ♦  &  tfi uk' *}-+&(rrg)«  i «•  $f) »*s, 

(4.2.1) 

with  Sd>  including  all  the  remaining  terms,  the  following  discussions  on 
discretization  will  be  based  upon  equation  (4.2.1). 


4.3  Discretization 


The  first  step  toward  discretization  is  the  arrangement  of  a  con¬ 
venient  grid  system.  For  the  present  study,  a  staggered  grid  was  adopted 
as  shown  in  Figures  4.1  and  4.2.  This  choice  minimizes  the  amount  of 
interpolation  needed  and  improves  stability.  The  locations  of  the  Rey¬ 
nolds  shear  stresses  were  also  staggered  to  the  places  as  indicated  in 
these  figures. 

Following  the  standard  approach  of  TEACH- type  codes,  equation 
(4.2.1)  is  integrated  over  a  control  volume  chosen  for  the  variable  <j>. 
Using  central  differencing  for  the  diffusion  terms,  the  integrated  form 
of  equation  (4.2.1),  with  all  the  subscripts  and  dimensions  referred  to 
Figure  4.3,  becomes 


“  Cwfar  *  ®  CwCej  *  Qg  ”$p) 

-D*  C 4>p-$w)  *  Dn  -$p)  -Ds  (#p  -$s)4 

where 

Ce  (fUx)e 


0,  -(raraB)*  (f  Ux)„ 


C*  -(razaS)* (fUt)n 


Cs  -(fAIA&)s  tfUr)s 


Q  -(Ar4z)d  (fUi)d 
Ci  -(Ar^z)*  (f d#)M 

Cfe  -lrarad)trt  /(Sz)t 
Dw  -(rAr46)Bn, /(tz)« 
D„  -(razAej.rWtSV), 
Ps  -(rAm6)tr»/(rr)j 
AV  "(rArAfiAz) 

S*  ”  av  js*  rdrdedz 


The  source  term  5$  A V  may  be  conveniently  separated  into  two  parts 
as  Sp$p  +  Sjj  in  order  to  take  full  advantage  of  the  solution  algorithm. 

A  more  elaborate  discussion  of  this  practice  and  its  benefits  will  be 
found  in  Patankar  [1980]. 

For  the  convection  terms,  a  proper  approximation  of  the  <|>'s  at  the 
control -volume  surfaces  has  been  proven  to  be  crucial  to  the  stability 
of  the  solution  algorithm.  The  accuracy  of  the  calculated  results  also 
depends  strongly  on  the  way  these  terms  are  simulated.  A  brief  outline 
is  given  next  of  the  schemes  chosen  for  the  present  study. 

1  For  <{>d  and  <j>u  which  stem  from  convection  along  the  main  flow  direc¬ 
tion,  upwind  differencing  was  invariably  used  with  *  <t>  and  <b  <j>°  where  $ 

w  H  ”  r 

stands  for  the  value  of  at  the  adjacent  upstream  station.  This  choice 
is  in  line  with  the  notion  of  a  parti  ally-parabolic  flow  in  that,  along 
this  direction,  there  is  no  flow  reversal  and  that  diffusion  is  negli¬ 
gible. 

2  For  <j>e,  $w,  $n  and  <j>s  over  a  cross-stream  plane,  two  types  of 
approximations  were  used:  the  central /upwind  (HYBRID)  scheme  and  the 
quadratic  upwind-weighted  interpolation  (QUICK)  scheme. 

HYBRID:  This  scheme,  first  proposed  by  Spalding  [1972],  Is  based  on 
an  approximation  to  the  exact  solution  of  the  one-dimensional  convec¬ 
tion-diffusion  equation 

4(eiV»-r$)«o 

between  two  points  with  known  4>'s.  It  employs  upwind  differencing  where 


the  flow  is  convection-dominated  and  shifts  to  central  differencing 
otherwise.  The  determining  parameter  is  the  local  Peclet  number  defined 
as  Pe  =  pU^  A£/T;  the  flow  is  considered  to  be  convection-dominated  if 
| Pe|  >  2.  Thus,  in  terms  of  an  approximation  for  <j>  ,  the  HYBRID  scheme 
yields 

4e  =  I  ($E  +  <PP)  if  |Pel  i  2 

4>e  *  <f>p  if  Pe  >  2 

4>e  4>e  if  Pe  <  -2. 

Similar  expressions  can  easily  be  constructed  for  the  <t>'s  at  the  other 
control -volume  surfaces.  Note  that  when  upwind  differencing  is  used 
at  the  condition  of  convection-dominance,  diffusion  at  the  corresponding 
cell  surface  is  neglected. 

It  has  long  been  known  that  upwind  schemes,  while  removing  the 
stability  problems  associated  with  central  differencing  schemes,  may 
introduce  severe  numerical  diffusion  arising  primarily  from  streamline- 
to-grid  skewness.  This  unphysical  diffusion  can  cause  serious  difficul¬ 
ties  for  turbulence  modelers  who,  in  order  to  evaluate  the  performance 
of  turbulence  models,  need  to  eliminate  from  their  computations  all  the 
inaccuracies  and  uncertainties  due  to  the  finite  difference  practices. 
Alternatives  to  the  HYBRID  scheme  have  been  proposed  and  tested  in  the 
literature,  see  Huang,  Launder  and  Leschziner  [1983]  for  a  brief  review. 

QUICK:  In  an  attempt  to  overcome  the  numerical  diffusion  problems 
associated  with  upwind  schemes  and  to  combine  the  accuracy  of  central 
differencing  with  the  stability  of  upwind  differencing,  Leonard  [1979] 
proposed  a  three-point,  upwind-weighted  quadratic  interpolation  to 
approximate  the  <j>'s  at  the  control -volume  surfaces.  The  proper  formulae 
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for,  say,  are 


4>e  =  +  +  ae^P  +  ^e 

(6^(1  +  5  - / o  ^  ^ )/8A. 

if 

ce 

>  0 

where 

a 

1-  V8*i-1 

if 

Ce 

<  0 

and 

/-  <«,)*<  Vsi + 

if 

Ce 

>  0 

if 

Ce 

<  0 

with  symbols  and  dimensions  defined  in  Figure  4.4.  Formulae  appropriate 
for  <j»'s  at  other  surfaces  are  summarized  in  Table  4.1.  Details  of  the 
QUICK  scheme  can  be  found  in  Leonard  [1979]  and  of  its  performance  in 
Han  et  al.  [1981]  and  Huang  et  al.  [1983]. 

With  the  formulae  given  above  and  the  difference  form  of  the  con¬ 
tinuity  relation 


equation  (4.3.1)  may  be  written  as 

ap4>p  s  3^4^  +  ®s^S  +  +  ®u^P  +  (4.3.2) 

where 

aP  *  aN  +  as  +  aE  +  aw  +  aU  "  SP  "  QP 
and  the  a's  denote  the  combined  contributions  of  convection  and  diffu¬ 
sion  at  the  control  volume  surfaces  to  the  balance  of  <pp  at  the  cell 
center. 

Depending  upon  the  differencing  scheme  chosen  for  the  <j>'s  over  a 
cross-stream  plane,  the  expressions  for  a's,  Qjj  and  Qp  take  the  different 


* 


forms  given  below. 

For  the  HYBRID  scheme: 


a£  =  ®”^e  *  ^e  ”  2  ’ 

aW  =  ®Cw  ’  Dw  +  T  »  °* 


aN  =  &“Cn  *  Dn  *  2°  * 


as  *  &cs  ’  °s  +  T  ’  01 


aU  =  Cu 


Qp =  Qu  =  0 


where  (  ]  stands  for  "the  largest  quantities  compared". 
For  the  QUICK  scheme: 


1 


aN  *  Dn  "  2  Cn 


aS  =  Ds  +  2  Cs 


aE  =  De  -  I  Ce 


aU  =  Dw  +  2  Cw 


aU  =  Cu 


Q_  =  C  a  -  C  a  +  C  a  *  C  a 
w  w  e  e  s  s  n  n 


Q. .  =  C  B  -  C  B  +  CB  -  CS 
wBw  e  e  s  s  nsn 


Equation  (4.3.2)  represents  the  general  form  of  the  finite  differ 


ence  equation  to  be  solved.  The  solution  procedure  adopted  for  the 
present  study  will  be  outlined  after  a  brief  description  in  the  next 
section  of  the  boundary  conditions  prescribed  numerically. 
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4.4  Numerical  Prescriptions  of  Boundary  Conditions 


Symmetry  of  the  flow  configuration  with  respect  to  the  plane  2=0 
allowed  the  performance  of  calculations  over  one  symmetrical  half  of  the 
duct  cross-section.  The  specification  of  required  boundary  conditions 
is  summarized  below. 


At  the  inlet  plane,  which  is  chosen  at  5  DH  upstream  of  the 


bend 


~7  ~7  _ 

section,  U  ,  U  ,  u  ,  u  and  u  u  were  interpolated  from  measurements 

*  y  *  y _  *  y 


at  the  same  station;  Uz>  uz  and  uxuz  were  deduced  from  symmetry  consider¬ 


ations;  and,  in  the  absnece  of  any  reliable  information,  uxuz  was  set  to 
zero.  The  distribution  of  c  was  prescribed  by  setting 

£  -  C  3/4  k3'2/* 
u  m 

with  i,m,  the  mixing  length,  found  from  a  generalization  of  the  straight 
pipe  formula  of  Nikuradse  [1932]  to  a  square  duct  geometry.  The  formula 
used  was 


£m  =  <Dh[Bz(1  -  K2  8z)  *  V1  -  1‘2  y]V2 

with 


S  »  I  (i  .  I] 
Z  D,,  u  D,.' 


and 


=  X 
dh 


0 


UH 


At  the  exit  plane,  chosen  at  5  DH  downstream  of  the  bend  exit, 
3Ux/3x  =  0  was  imposed. 

Along  the  symmetry  plane  at  z  «  0,  the  derivatives  of  U  ,  U  ,  P, 

x  y 

k  and  t  as  well  as  Uz  itself  were  set  to  zero. 

At  solid  walls,  the  wall-function  approach  outlined  in  Launder  and 


Spalding  [1974]  was  adopted.  This,  in  essence,  consists  in  specifying 
boundary  conditions  near  the  walls  rather  than  at  the  walls.  The  first 


calculation  node,  P,  is  placed  at  a  distance  yp  away  from  a  wall  in  a 
region  where  the  logarithmic  law  of  the  wall  prevails  and  in  which  the 
turbulence  is  in  a  state  of  local  equilibrium  and  the  shear  stress  r 

P 

is  approximately  equal  to  that  of  the  wall,  x^.  Within  this  equilibrium 
layer,  it  can  be  shown  that 

t„  -*p<^%**Up/AaEyf  \4U/9) 

kp  *  Ut  / And  Cp  ■  UiVxyp 

wkcre 

Uc  -  lx./p)1/2  4nd  E  -  9.793  . 

These  were  the  conditions  imposed  near  wall.  Since  there  are  no  relia¬ 
ble  experimental  correlations  for  the  form  of  the  law  of  the  wall,  if 
exists  at  all,  for  complex,  three-dimensional  skewed  boundary  layer 
flows,  the  lateral  velocities  Ur  and  Uz  were  assumed  to  follow  the  same 
logarithmic  relation  as  that  for  IL  outlined  above. 

4.5  Solution  Algorithm 

In  the  present  study  the  elliptic  numerical  procedure  of  Humphrey 
[1977]  was  modified,  according  to  the  guidelines  of  Pratap  [1975],  to 
perform  parti  ally-parabolic  calculations.  The  pertinent  features  of 
the  solution  algorithm  are: 

1.  Calculation  is  performed  by  marching  through  the  flow  domain  along 


the  main  flow  direction  as  many  times  as  is  required  for  a  predeter 
mined  convergence  criterion  to  be  achieved. 

2.  The  pressure  distribution  is  always  stored  in  a  three-dimensional 
array  for  the  entire  flow  field  so  that  it  can  be  corrected  and 
updated  during  each  sweep  of  the  flow  domain. 

3.  The  remaining  flow  variables  are  continuously  recalculated  and  are 
stored  temporarily  in  two-dimensional  arrays  at  the  current  com¬ 
puting  station  and,  for  velocity  components  only,  at  the  correspond 
ing  upstream  station.  Because  of  this  practice,  nonlinear  convec¬ 
tion  terms  in  the  momentum  equations  are  linearized  with  respect 

to  their  values  at  the  upstream  station  (rather  than  their  values 
at  the  current  station  from  the  previous  sweep). 

4.  The  hydrodynamic  variables  are  solved  by  a  slightly  modified  ver¬ 
sion  of  the  SIMPLE  procedure.  For  a  detailed  account  of  the  SIMPLE 
procedure,  see  Patankar  [1980], 

5.  The  discretized  equations  are  solved  by  a  line-by-line  iterative 
procedure— the  tridiagonal  matrix  algorithm  (TDMA) . 

The  major  calculation  steps  are  summarized  below: 

1.  Assign  initial  guessed  values  to  the  pressure  field. 

2.  Solve  the  momentum  equations  for  the  cross-stream  velocity  compon¬ 
ents  U  and  II  . 

r  z 

3.  Solve  the  Ug  momentum  equation.  Note  that  UQ  is  located  ahead  of 
the  plane  containing  Ur  and  Uz  (see  Figure  4.2). 

4.  Update  the  pressure  and  velocity  distributions  according  to  the 
modified  SIMPLE  algorithm. 

5.  Solve  the  remaining  transport  equations  (k  and  e). 


6.  Calculate  the  Reynolds  stresses  from  the  ASM  relations  (ASM  only). 

7.  March  to  the  next  downstream  station  and  repeat  Steps  2  to  6  until 
the  exit  plane  of  the  calculation  domain  is  reached. 

8.  Repeat  Steps  2  to  7  with  the  most  recent,  more  correct  pressure 
distribution  as  the  new  initial  guess  until  convergence  is  achieved. 
The  immediate  advantage  of  the  partially-parabolic  procedure  over 

the  elliptic  one  is  that  the  storage  requirement  for  the  same  grid  den¬ 
sity  is  drastically  reduced  for  the  former.  Due  to  this  tremendous 
savings  on  computer  memory,  it  becomes,  in  principle,  feasible  to  compute 
three-dimensional  flows  on  a  grid  system  refined  enough  to  render  numeri¬ 
cal  diffusion  an  insignificant  portion  of  the  numerical  accuracy.  For  a 
more  detailed  exposition  of  the  partially-parabolic  procedure,  see  the 
Ph.D.  dissertation  of  Pratap  [1975]. 


CHAPTER  5.  VALIDATION 


5.1  Introduction 


The  finite-difference  forms  of  the  governing  equations,  together 
with  the  appropriate  boundary  conditions  were  coded,  following  the 
guidelines  of  Pratap  [1975],  into  a  computer  program  which  solves  three 
dimensional,  partially-parabolic  flows  in  cylindrical/rectangular  coor¬ 
dinates.  Prior  to  predicting  the  main  flow  of  interest  to  the  present 
study,  rigorous  testing  was  performed  against  several  flows  with  either 
known  analytical  or  reliable  numerical  solutions  or  generationally 
accepted  experimental  information. 

The  purpose  of  testing  was  twofold: 

1 .  To  check  the  correctness  of  the  procedure. 

The  present  algorithm  was  obtained  by  modifying  the  eliiptic  pro¬ 
cedure  of  Humphrey  [1977].  In  order  to  ensure  that  the  modifications 
had  been  carried  out  properly,  predictions  of  laminar  flows  in  the  fol¬ 
lowing  geometries  were  performed: 

i  two-dimensional  straight  channel; 
ii  two-dimensional  curved  channel; 

iii  straight  duct  of  square  cross-section  with  one  wall  moving  at 
a  constant  speed  normal  to  the  main  flow  direction; 
iv  90°  bend  of  square  cross-section  with  a  long  straight  duct 
preceeding  it. 


The  results  of  the  testing  are  summarized  and  discussed  in  Section 


2.  To  evaluate  the  performance  of  the  k-c  model  and  the  ASM 
closures  adopted  for  this  study. 

Numerous  examples  of  predictions  obtained  with  these  two  models 
have  already  been  published  in  the  literature  (see,  for  example,  the 
review  of  Rodi  [1978]).  The  present  tests  were  conducted  to  ensure 
that  there  were  no  coding  errors  in  the  portion  of  the  program  concern¬ 
ing  the  turbulence  models  and,  more  importantly,  to  evaluate  the  per¬ 
formance  of  these  models  within  the  framework  of  a  partially-parabolic 
procedure.  The  test  geometries  selected  were  a  two-dimensional  straight 
channel  and  a  two-dimensional  curved  channel.  The  results  of  these  tur¬ 
bulent  flow  tests  will  be  discussed  in  Section  5.3. 

5.2  Laminar  Flow  Tests 

In  this  section,  the  four  laminar  flow  test  cases  are  described  and 
their  results  briefly  discussed.  The  two-dimensional  flow  tests  of  i 
and  ii  were  predicted  by  imposing  two  symmetry  conditions  along  the 
third  (z)  direction  of  the  three-dimensional  code.  Calculations  of 
cases  i,  ii  and  iii  were  carried  out  by  using  a  plug  flow  profile  as 
the  inlet  condition,  whereas  a  fully  developed  velocity  profile  was 
imposed  at  the  inlet  plane  in  case  iv. 

Case  i  Developing  flow  in  a  two-dimensional  straight  channel 

A  HYBRID  scheme  calculation  on  a  (y  =  17)  x  (x  =  101)  equally 
spaced  grid  was  performed  over  a  symmetrical  half  of  the  channel.  As 
shown  in  Figure  5.1  the  calculated  fully  developed  velocity  profile 
follows  the  analytical  solution  almost  identically,  with  the  maximum 


deviation,  which  occurs  at  the  first  node  next  to  the  channel  wall, 
being  less  than  2%. 

Case  ii  Developing  flow  in  a  two-dimensional  strongly  curved  channel. 

In  order  to  test  the  appropriateness  of  the  concept  of  parti  ally- 
parabolic  flows,  a  strongly  curved  channel  was  purposefully  chosen. 
Figure  5.2  shews  the  calculated  fully  developed  velocity  distribution 
as  compared  to  the  exact  solution  (Goldstein  [1965]).  Agreement  is 
excellent  and  the  maximum  deviation  from  the  analytical  solution,  which 
occurs  at  the  first  node  next  to  the  inner  wall  at  r  =  r. ,  is  less  than 
3.5%.  This  calculation  was  performed  with  the  HYBRID  scheme  using  a 
(r  =  22)  x  (9  =  183)  uniformly  spaced  grid. 

Case  iii  Developing  flow  in  a  straight  duct  of  square  cross-section 
with  one  wall  sliding  normal  to  the  main  flow  direction  at  a  constant 
speed. 

This  flow  was  calculated  with  both  the  HYBRID  and  the  QUICK  schemes 
on  an  equally  spaced  grid  of  (y  *  15)  x  (z  =  15)  x  (x  =  121).  The  cal¬ 
culation  of  Burggraf  [1966]  was  chosen  as  the  standard  for  comparison. 
Figure  5.3  demonstrates  clearly  the  superiority  of  the  QUICK  scheme  over 
the  HYBRID  scheme.  In  fact,  a  calculation  with  the  QUICK  scheme  using  a 
(y  *  8)  x  (z  =  8)  x  (x  *  121)  grid  produced  results  similar  to  those 
generated  in  the  aforementioned  HYBRID  calculation. 

Case  iv  Flow  in  a  90°  bend  of  square  cross-section  with  fully  developed 
velocity  profile  close  to  the  inlet  of  the  bend. 

The  particular  flow  calculated  corresponds  to  the  measurements  of 
Humphrey  et  al .  [1981]  at  Re  *  790  and  De  *  368.  Both  the  HYBRID  and 


the  QUICK  schemes  were  used  on  a  (r  =  17)  x  (z  =  11 )  x  (6  =  36)  uniformly 
spaced  grid.  The  agreement  with  measurements  was  generally  very  good  at 
bend  angles  of  0°,  30°  and  90°,  with  discrepancies  of  U„  up  to  30%  being 
found  at  9  =  60°.  Figure  5.4  compares  the  calculated  longitudinal  velo¬ 
city  profiles  at  Z  =  0.5  and  9  =  90°  (refer  to  Figure  2.1  for  coor¬ 
dinate  system)  with  experimental  results.  The  QUICK  scheme  calculation 
again  yielded  better  overall  agreement  with  the  measured  data. 

5.3  Testing  of  Turbulence  Models 

The  laminar  flow  tests  presented  in  the  previous  section  have 
demonstrated  the  applicability  of  the  parti ally-parabol ic  procedure 
for  predicting  internal  flows  of  characteristics  similar  to  the  present 
study.  This  section  provides  further  validation  on  the  turbulence 
models  used. 

Case  i  Flow  in  a  two-dimensional  straight  channel. 

This  is  the  simplest  possible,  yet  very  informative  test  that  can 
be  performed  by  using  the  current  computer  program  with  minor  modifica¬ 
tions.  For  the  purpose  of  comparison,  the  experiment  of  Laufer  [1950] 
at  Re  *  123,200  was  simulated  on  an  unequally  spaced  grid  of  (y  =  18)  x 
(x  =  241)  over  a  symmetrical  half  of  the  channel.  Following  the  experi¬ 
mental  correlation  of  Laufer,  the  log-wall  constants  of  <  =  0.33  and 
E  =  6.14  were  used  instead. 

Figures  5.5  and  5.6  show  the  calculated  distributions,  using  k-e 
model,  of  the  mean  longitudinal  velocity,  Ux>  and  the  turbulence  energy, 
k,  as  compared  with  the  experimental  measurements.  A  preliminary  check 


on  the  "correctness"  of  the  ASM  was  made  by  substituting  the  calculated 
values  of  Ux,  k  and  e  into  the  ASM  relations  given  in  Section  3.3.3  and 
solving  for  the  Reynolds  stresses.  It  can  be  seen  from  Figure  5.7  that, 
assuming  U^,  k  and  e  as  calculated  by  the  k-e  model  were  correct,  the 
ASM  closure  is  capable  of  predicting  broadly  satisfactory  distributions 
of  the  Reynolds  stresses. 

Direct  computation  using  the  ASM  closure,  however,  yielded  mean 
velocity  and  Reynolds  stresses  in  poor  agreement  with  the  measurements. 

In  particular,  the  level  of  u^uj  was  twice  as  high  compared  to  the  data. 
The  alternative  approach  of  employing  the  isotropic  k  and  e  equations 
[equations  (3.3.1.11)  and  (3.3.1.12)],  as  has  been  done  by  many  of  the 
ASM  users,  produced  better  results.  Figures  5.8  to  5.10  show  that  the 
general  trend  of  the  predictions  were  satisfactory,  with  the  most  serious 
discrepancies  being: 

1.  Very  close  to  the  wall  the  value  of  ITiT’  falls  off  instead  of 
rising  as  suggested  by  the  data. 

2.  The  differences  between  ulT  and  TTu^"  are  too  large. 

In  order  to  check  the  sensitivity  of  the  model  to  the  model  con- 

I  I 

stants,  a  further  run  with  *  0.35  and  Z^  -  0.20  was  carried  out. 
Figures  5.8  to  5.10  show  that  better  agreement  was  achieved  with  the  new 
constants,  especially  for  the  distributions  of  u  u  and  U  .  Since  the 

it  y  y  x 

same  change  of  C-j  and  C2  did  not  produce  appreciable  overall  improvement 

I 

in  the  2D  curved  channel  case,  the  originally  recommended  values  of  C1  = 

I 

0.5  and  C2  =  0.3  were  used  in  the  final  computations  described  in 
Chapter  6. 


Case  ii  Flow  in  a  two-dimensional  curved  channel. 

The  development  of  turbulent  shear  flows  on  curved  surfaces  is 
characterized  by  a  strong  sensitivity  of  the  turbulence  structure  to 
streamline  curvature  (Bradshaw  [1973]).  Attempts  to  predict  such 
flows  by  using  eddy  viscosity  transport  closures,  such  as  the  k-e 
model,  have  not  been  very  successful  unless  empirical  means  to  account 
for  curvature  effects  are  incorporated  into  the  models.  Approaches  of 
this  type  have  been  reported  by,  among  others,  Launder  et  al.  [1977] 
and  Humphrey  and  Pourahmadi  [1982]  with  some  success  for  moderately 
curved  channels.  It  is  argued  (Gibson  [1978])  that  only  turbulence 
models  based  upon  the  calculation  of  Reynolds  stresses  directly  from 
their  transport  equations  can  account  for  the  streamline  curvature 
effects  correctly. 

In  order  to  test  the  capability  of  the  ASM  closure  for  predicting 

the  curvature  effects  mentioned  above,  the  measurements  of  Eskinazi 

and  Yeh  [1956]  in  a  curved  rectangular  duct  of  large  aspect  ratio  was 

simulated  on  a  (r  *  24)  x  (9  =  184)  unequally  spaced  grid.  Figures 

5.11  and  5.12  compare  the  calculated  (HYBRID  scheme  only)  distributions 

of  U0,  UgUg,  u^  and  u^g  with  the  experimental  data.  While  the  ASM 

closure  tends  to  overpredict  the  turbulence  intensities  uTuT  and  u  u 

oBrr 

near  the  convex  (stabilizing)  wall  and  underpredict  the  values  of  UQ  and 
urUg  near  the  concave  (destabilizing)  wall,  predictions  are,  in  general, 
in  good  agreement  with  the  data  of  Eskinazi  and  Yeh. 

The  idea  of  using  a  non-symmetric  wall -proximity  function  f  [equa¬ 
tions  (3.3.2.6)-(3.3.2.8)],  first  proposed  by  Humphrey  and  Pourahmadi 
[1982],  was  also  tested.  This  correction  produced  very  little  changes 


on  the  distribution  of  IL/IL  and  its  effect  on  uQuQ  and  u  u  is  shown 

o  u ,max  u  y  r  r 

in  Figures  5.11  and  5.12.  Since  this  non-symmetric  f-function  did  not 
yield  results  in  better  overall  agreement  with  the  data,  the  original 
(symmetric)  f-function  [equation  (3. 3. 2. 8)]  was  retained. 


CHAPTER  6.  RESULTS  AND  DISCUSSION 


6.1  Introduction 

In  this  chapter,  detailed  measurements  of  developing  turbulent  flow 
in  a  square  cross-section  180°  bend  and  its  upstream  and  downstream 
tangents  are  presented  and  discussed.  These  data  are  compared  with 
corresponding  predictions  by  both  the  k-e  model  and  the  algebraic  stress 
model  at  selected  streamwise  stations.  Section  6.2  considers  the 
measurements  and  the  k-e  model  predictions.  Both  the  HYBRID  and  the 
QUICK  schemes  were  used  in  the  calculations.  When  QUICK  and  HYBRID 
results  coincide,  a  single  profile  is  shown  in  the  figures.  Otherwise, 
HYBRID  results  are  plotted  as  continuous  lines  and  QUICK  as  dashed 
lines.  It  should  be  noted  that  dashed  lines  have  also  been  used  in 
some  graphs  for  plotting  best  fits  to  experimental  data  with  a  larger 
than  normal  degree  of  scatter  at  the  2z/DH  =0.75  and  0.875  spanwise 
locations.  Section  6.3  presents  the  predictions  with  the  ASM  closure. 
Since  it  has  not  been  successful  in  producing  converged  results  with 
QUICK  scheme  at  the  time  of  writing,  only  HYBRID  calculations  are  in¬ 
cluded.  Two  simplified  "experimental"  tests  of  the  ASM  closure  are 
also  presented  in  this  section. 

6.2  Turbulent  Flow  Measurements  and  the  k-e  Model  Calculations 

The  calculation  presented  in  this  section  was  performed  on  an 
unequally  spaced  grid  consisting  of  the  following  distributions  of  nodes 


(z  =  14)  x  (r  =  25)  x  (e  =  45)  In  the  bend,  and  (z  =  14)  x  (r  =  25)  x 
(x  =  20)  in  the  upstream  and  downstream  tangents,  both  of  length  5  D^. 

The  grid  distribution  over  a  cross-stream  plane  is  listed  in  Table  6.1. 

The  convergence  criterion  was  that  the  maximum  of  the  normalized  resi- 
dual  summations  at  every  cross-sectional  plane  be  less  than  10  . 

Measurements  of  the  flow  in  the  upstream  tangent,  taken  at  XH  =  -5 
and  -1  are  shown  in  Figures  6.3  and  6.4  respectively.  Comparing  the 
data  at  these  two  stations,  particularly  the  Reynolds  stress  measure¬ 
ments  shown  in  Figure  6.14,  suggests  that  the  flow  is  still  developing 
at  =  -5,  after  travelling  a  distance  of  30  from  the  uni formi zing 
section.  Figures  6.2  and  6.14  provide  comparisons  of  the  measurements 
at  XH  =  -1  with  corresponding  profiles  interpolated  from  the  data  ob¬ 
tained  by  Melling  and  Whitelaw  [1976]  at  36.8  hydraulic  diameters  in  a 
straight  duct  of  square  cross-section.  The  two  data  sets  are  in  agree¬ 
ment  to  within  the  experimental  error  of  the  measurements.  Since  the 
measurements  by  Melling  and  Whitelaw  correspond  to  an  essentially  devel¬ 
oped  flow,  the  differences  shown  in  Figure  6.2  for  the  radial  (transverse) 
velocity  component  are  attributed  mainly  to  the  elliptic  influence  of 
the  bend  on  the  flow  in  the  tangent.  Similar  observations  on  the  in¬ 
fluence  of  a  bend  on  its  upstream  tangent  flow  have  been  noted  by  Hum¬ 
phrey  et  al.  [1981]  and  Taylor  et  al.  [1982]  in  the  same  90°  bend  con¬ 
figuration  but  with  differing  Inlet  conditions.  Taylor  et  al.  found 
that  in  a  turbulent  flow  with  relatively  thin  boundary  layers  the  bend 
influenced  the  measurements  taken  at  XH  =  -0.75,  their  furthest  upstream 
location..  Measurements  at  XH  =  -2.5  taken  by  Humphrey  et  al.  with  essen¬ 
tially  fully  developed  flow  conditions  in  the  upstream  tangent  suggested 


only  a  weak  elliptic  effect  of  the  bend  on  the  flow  at  this  location. 
However,  these  authors  did  not  measure  the  transverse  velocity  compon¬ 
ent;  their  commentary  is  based  on  the  unchanged  appearance  of  stream- 
wise  velocity  and  turbulence  intensity  contours  between  XH  =  -1  and 
=  -2.5.  By  contrast,  for  the  conditions  of  the  present  experiment, 
the  velocity  data  suggests  that  the  favorable  streamwise  pressure 
gradient  along  the  inner  wall  of  the  bend  (see  Figure  6.1)  induces  a 
mean  transverse  flow  directed  at  theinnerwall  which  is  already  notice¬ 
able  at  X^  «  -5.  Although  weak  (Ur/Ug  =  0.02),  it  appears  that  the 
transverse  flow  induced  by  pressure  forces  in  the  upstream  duct  over¬ 
comes  the  weaker  stress-driven  cross-stream  motion  (shown  clearly  in 
the  measurements  of  Melling  and  Whitelaw),  and  results  in  Ur  profiles 
whose  shapes,  and  variations  of  shape  with  z  location,  agree  qualita¬ 
tively  with  a  simple  superposition  of  the  pressure-induced  and  stress- 
driven  cross-stream  flows. 

Calculations  of  Uf  at  e  *  3°  (and  between  this  location  and  X^  =  -5, 
not  given  here)  always  showed  the  flow  moving  towards  the  convex  wall- 
side  of  the  test  section.  Because  the  turbulence  model  is  insensitive 
to  stress-driven  secondary  motions,  this  result,  and  the  calculated 
displacement  of  the  maximum  in  UQ  towards  the  convex  wall,  clearly 
support  the  experimental  finding  that  elliptic  effects  are  transmitted 
from  the  bend  into  the  upstream  tangent  via  the  pressure  field,  farther 
than  has  previously  been  ovserved.  As  with  the  Ur  velocity  component, 
the  discrepancies  shown  between  calculated  and  measured  turbulent 
stresses  at  9  *  3°  arise  from  the  assumption  of  Isotropy  in  the  tur¬ 
bulence  model.  To  predict  more  accurately  the  cross-stream  motion  and 


turbulent  stresses  in  the  upstream  tangent  requires  modeling  of  important 
anisotropic  flow  characteristics  such  as  near-wall  pressure-strain 
effects.  These  effects  cannot  be  dealt  with  by  isotropic  viscosity 
models  of  any  kind. 

The  mean  velocity  and  turbulent  stress  data  taken  at  different 
streamwise  locations  in  the  bend  are  plotted  in  Figures  6. 5-6. 9  and 
6.15.  The  pressure  drop  measured  at  the  side  wall  of  the  bend  is  shown 
in  Fig.  6.1.  In  general,  the  sense  of  the  flow  up  to  0  =  90°  is  in 
agreement  with  earlier  observations  of,  for  example,  Humphrey  et  al. 
[1981]  and  Taylor  et  al.  [1982].  Closer  inspection  of  the  plots  reveals 
additional  interesting  features.  The  measurements  of  the  pressure  coef¬ 
ficient  Cp  show  the  opposing  pressure  gradient  initially  expected  at  the 
concave  wall  and  the  favorable  gradient  at  the  convex.  In  contrast  to 
the  data  measured  by  Taylor  et  al.  for  a  90°  bend,  the  maximum  and 
minimum  values  of  Cp  do  not  coincide  at  the  same  streamwise  location. 

In  this  work  the  value  of  Cp  maximizes  at  the  concave  wall  at  0  =  45°, 
and  attains  its  minimum  value  at  the  convex  wall  at  0  =  177°.  As  of 
0  *  45°,  the  streamwise  pressure  gradient  is  favorable  throughout  most 
of  the  bend.  At  0  *  3°  the  streamwise  velocity  profiles  all  shown  their 
maximum  values  displaced  towards  the  inner  radius  wall,  due  to  the 
favorable  streanwise  pressure  gradient  there.  The  radial  component  of 
velocity  is  everywhere  directed  towards  the  inner  wall  except  for  a 
small  flow  region  about  the  symmetry  plane.  This  region  marks  the  incep¬ 
tion  of  secondary  flow,  driven  by  the  transverse  pressure  gradient  which 
arises  due  to  lateral  curvature  of  the  main  flow  in  the  bend.  The  secon¬ 
dary  flow  is  more  intense  at  0  =  45®  as  shown  in  Figure  6.6.  At  0  =  45°, 


the  profile  for  Ur  near  the  side  wall  (2z/DH  =  0.75)  shows  a  large  nega¬ 
tive  velocity,  while  the  profiles  at  the  other  z  locations  are  large  and 
positive.  The  sense  of  motion  is  from  the  convex  to  the  concave  surface 
along  the  symmetry  plane,  and  back  to  the  convex  surface  along  the  side 

walls  of  the  bend.  As  of  6  =  45°,  measurements  of  U  at  all  z  locations 

r 

show  this  component  always  positive  (directed  from  the  convex  to  the 

concave  surface).  This  means  that  the  cross-stream  plane  return  flow 

adjacent  to  a  flat  wall  in  the  bend  is  confined  to  a  narrow  region  less 

than  D^/8  in  width.  A  simplified  mass  balance  at  6  *  90°  suggests  that 

the  radial  component  of  velocity  in  this  narrow  region  varies  between  0 

at  the  wall  and  a  maximum  of  about  0.40  x  Ug  at  the  peak  location.  This 

result  is  in  qualitative  agreement  with  the  observations  of  Taylor  et 

al.  [1982].  The  highest  return  flow  measured  by  them  (0.40  x  Ug)  was 

found  at  2z/DH  =0.95  and  8  =  60°  in  a  90°  bend  with  Re  =  40,000  and 

Rc/Dh  =  2.3.  Calculations  of  the  present  flow  gave  a  peak  value  of 

Ur/UB  =  0.30  at  2z/Dh  =  0.95,  at  the  90°  plane. 

Between  0=3°  and  90°,  the  influence  of  destabilizing  curvature 

raises  the  levels  of  all  the  measured  Reynolds  stresses  at  the  concave 

wall  of  the  bend.  The  effect  is  particularly  noticeable  in  the  plots 

for  the  stresses  at  6  =  45°,  and  has  decreased  by  the  time  the  flow 

reaches  the  0  *  90°  plane.  Similar  observations  have  been  made  by 

Eskinazi  and  Yeh  [1956]  in  a  two-dimensional  curved  channel  flow,  and 

by  Humpnrey  et  al.  [1981]  and  Taylor  et  al.  [1982]  in  their  respective 

curved  duct  flows.  In  particular,  Eskinazi  and  Yeh  demonstrate  the 

“2 

strong  generation  of  ur  near  a  concave  wall  and  its  correspond  suppression 
near  a  convex  wall.  While  at  0  *  45°  the  present  flow  is  already  three- 


dimensional,  the  results  for  ur/Ug  are  in  qualitative  agreement  with  the 
observations  of  Eskinazi  and  Yeh.  Destabilizing  and  stabilizing  curva¬ 
ture  effects  at  the  respective  concave  and  convex  walls  of  the  present 

flow  are  responsible  for  producing  large  levels  of  anisotropy.  For 

1 — T~ 

example,  at  0  =  45°,  u  „/u  „  =  2  near  the  concave  wall.  Wall -dampening 

u  r 

of  radial  fluctuations  in  the  flow,  coupled  with  pressure  redistribution 

and  turbulent  diffusion  of  energy  between  the  normal  stresses,  account 

~2~ 

for  the  larger  levels  of  u  .  . 

y 

At  9  =  90°  and  130°,  plots  of  the  two  velocity  components  and  their 
corresponding  normal  stresses  show  striking  variations  in  the  radial 
coordinate  direction.  These  take  the  form  of  relatively  large  decelera¬ 
tions  in  the  mean  flow  components  at  about  n  *  0.4  and  are  accompanied 

by  relatively  large  increases  in  tL  and  decreases  in  iL  respectively  at 

y  r 

the  same  locations.  The  streamwise  deceleration  of  IL  and  U  is  due  to 

9  r 

the  "pumping"  of  low  speed  fluid  from  the  peripheral  region  of  the  duct 
into  the  core  of  the  flow.  Whereas  viscous  effects  provide  the  mechanism 
for  flow  retardation  at  the  walls,  it  is  the  inviscid  mechanism  of 
lateral  flow  curvature  which  drives  the  cross-stream  motion.  The  phen¬ 
omenon  has  been  analyzed  and  described  in  depth  by  Hawthorne  [1951],  and 
Horlock  and  Lakshminarayana  [1973]  discuss  it  extensively  in  the  context 
of  turbomachinery  aerodynamics.  In  particular,  the  latter  authors  pro¬ 
vide  general  expressions  for  the  generation  of  streamwise  vorticity 
including  the  influence  of  Bernouilli  surface  rotation  and  viscous 
effects.  Hawthorne  shows  that  In  flows  where  the  angle,  <J>,  between 
the  direction  of  normal  acceleration  and  the  normal  to  the  Bernoulli 
surface  (a  surface  of  constant  total  pressure)  is  not  0,  streamwise 


vorticity  is  induced.  For  a  small  bend  deflection  angle  ?  in  a  flow  with 


initial  uniform  vorticity  fi  in  the  plane  of  the  bend,  the  expression 
derived  by  Hawthorne  [1951]  is: 

5e  *  -  2  “o  6 

The  same  result  was  obtained  by  Squire  and  Winter  [1951)  using  a  differ¬ 
ent  analytical  approach  and  is  valid  for  Bernoulli  surfaces  which 
remain  undisturbed  during  passage  through  bend  (<j>  =  constant). 

It  is  important  to  realize  that  the  inviscid  mechanism  acts  over 
the  entire  cross-section  of  the  flow  (provided  $  t  0).  For  the  kind  of 
flow  of  interest  here,  the  distortion  of  Bernoulli  surfaces  during 
passage  of  the  flow  around  the  bend  cannot  be  neglected.  A  first 
approximation  accounting  for  this  distortion  leads  to  an  equation  for 
a,  the  angle  by  which  the  Bernoulli  surfaces  turn.  The  result  derived 
by  Hawthorne  [1951]  is: 

^-Sf  *  cosa 
de* 

This  equation  is  analogous  to  that  governing  the  motion  of  a  pendulum 
and,  as  shown  by  Hawthorne,  predicts  that  the  flow  through  a  bend  should 
oscillate  between  a  s  0  and  a  *  i?  with  a  period  (or  bend  angle)  for  a 
complete  oscillation  approximately  equal  to  2tt //  RC/DH'.  The  secondary 
motion  is  analogous  to  the  kinetic  energy  of  the  pendulum  and  also 
oscillates  with  the  same  period,  passing  through  zero  after  each  itvDh/Rc 
radians  of  turn.  For  the  present  flow  a  complete  oscillation  corres¬ 
ponds  to  about  197°.  This  means  that  at  about  98°  the  inviscid  mechanism 


starts  working  to  oppose  the  original  sense  of  the  cross-stream  motion, 

with  a  maximum  negative  amplitude  at  9  =  148°.  The  result,  as  shown  in 

Figures  6.8  and  6.9  for  e  =  130°  and  177°,  is  for  high  speed  fluid  to 

be  restored  to  the  core  of  the  flow. 

For  9  90°  and  130°,  plots  of  the  Reynolds  stresses  show  large 

changes  at  the  radial  positions  where  IL  and  U  have  been  decelerated. 

y  r 

Transport  equations  for  the  turbulent  stresses  can  be  obtained  (Bryant 

T  T 

and  Humphrey  [1976]),  and  for  uA  and  u  they  take  the  forms  given  below: 
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In  the  above  equations,  p!  and  represent  pressure  strain  redistribu- 

U  I 


tion  terms,  D.  and  D  denote  the  effects  of  turbulent  diffusion  and 
o  r 


dissipation,  and  viscous  diffusion  is  neglected.  The  terms  written  out 
explicitly  on  the  right-hand  side  of  these  two  equations  represent 
generation  of  the  stress  component  and,  hence,  of  the  kinetic  energy  of 
turbulence,  k.  Analysis  of  the  normal  stress  equations,  including  the 
effects  due  to  pressure  strain,  turbulent  diffusion  and  dissipation, 
while  desirable,  is  hampered  by  the  unavailability  of  appropriate  experi¬ 
mental  data;  hence,  the  contributions  of  the  latter  terms  to  the  respec¬ 
tive  balances  cannot  presently  be  established.  Nevertheless,  simplified 
analysis  of  the  generation  terms  alone  provides  a  basis  for  checking 
consistency  in  the  meausrements  and  can  shed  light  on  the  behavior  of 
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the  flow.  In  the  vicinity  of  the  bend  symmetry  plane,  symmetry  consider¬ 
ations  support  the  notion  that  the  generation  terms  in  equations  (6.2.1) 
and  (6.2.2)  involving  the  z  coordinate  direction  should  be  small  rela¬ 
tive  to  those  involving  variations  in  the  r  and  6  directions.  Approxi- 
~2  T 

mate  balances  for  uA  and  u  are  then  given  by: 


A  B  C 


D  E  F 


At  90°  uQu„,  is  positive  everywhere  and,  except  for  a  small  region 
w  r 

0.2  £  n  <  0.4  where  3U0/3r  is  negative,  term  A  represents  a  negative 

~2 

contribution  to  the  balance  of  uQ.  From  0  *  45°  to  90°  the  streamwise 

U 

velocity  component  U.  is  decelerated  between  n  =  0.2  and  0.5  by  the 

U 

inviscid  mechanism  explained  above.  The  result  is  for  term  B  in  equa¬ 
tion  (6.2.3)  to  contribute  positively  (together  with  term  A)  to  the 
~7 

balance  of  uQ.  At  all  radial  locations  near  the  symmetry  plane,  term  C 

T 

represents  a  reduction  in  the  magnitude  of  u0,  but  because  of  the  small 

values  of  the  ratio  U /r,  the  contribution  of  this  term  relative  to  A 

r 

and  B  is  small.  As  shown  in  Figure  6.7  ,  the  net  result  is  to  produce 
a  local  increase  in  the  magnitude  of  CL  between  n  8  0.4  and  0.5  approxi- 

U 

mately. 

In  a  similar  manner,  the  localized  minimum  in  u  at  n  -  0.4  in  the 

r 

90°  plane  can  be  explained.  At  all  radial  locations  term  E  in  equation 

— j r 

(6.2.4)  contributes  positively  to  the  balance  of  ur,  although  weakly  near 
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the  position  of  the  minimum  in  because  of  the  small  magnitude  of 

unu  there.  Between  0  =  45°  and  90°  the  plots  for  U  show  that,  as  for 

UQ,  U  has  also  been  decelerated  leading  to  a  positive  contribution  of 
6  r 

term  F  to  the  balance  of  uf.  However,  the  relatively  strong  streamwise 
deceleration  in  Ur  induces  large  radial  variations  in  its  profiles  pro¬ 
ducing  regions  in  the  flow  where  3Ur/3r  is  large  and  positive.  As  shown 
in  Fig.  6.7,  this  occurs  between  0  =  0.2  (at  2z/D^  =  0.75)  and  0.7  (at 
2z/Dh  =  0)  and  leads  to  a  negative  contribution  of  term  D  to  the  balance 

~2  ~ 
of  ur.  Since  term  D  is  the  only  one  contributing  negatively  to  it 

~~2 

must  be  the  cause  for  the  local  minima  in  the  profiles  for  ur  at  0  =  90°. 

Therefore,  the  effect  of  this  term  in  the  balance  must  be  large.  It 

appears  that  the  two  terms,  B  in  equation  (6.2.3)  and  D  in  equation 

~2  ~2 

(6.2.4),  are  the  principle  source  and  sink  for  uQ  and  ur  respectively 

in  this  simplified  analysis.  The  gradual  disappearance  of  the  maxima  in 

un  and  the  minima  in  Ci  between  9  =  130°  and  177°  is  due  to  the  inviscid 
u  r 

oscillatory  nature  of  the  flow  since,  by  restoring  high  speed  fluid  to 

the  core  of  the  flow  (see  Figures  6.8  and  6.9),  the  respective  roles  of 

terms  B  and  D  in  equations  (6.2.3)  and  (6.2.4)  are  reversed. 

"T  “7 

The  anistropic  effects  in  u.  and  u„  discussed  above  can  not  be 

y  r 

reproduced  by  the  calculations.  Nor  is  it  surprising  that  the  minima 
observed  in  the  experimental  profiles  for  the  mean  velocity  components 
at  9  ■  90°  and  130°  are  not  predicted.  The  assumption  of  local  isotropy 
in  the  turbulence  model  produces  such  large  levels  of  false  physical 
diffusion  as  to  preclude  an  accurate  spatial  resolution  of  the  flow. 

For  example,  at  9  *  90°  in  the  symmetry  plane  of  the  bend,  estimates  of 
v^/v  from  the  measurements  suggest  that  the  ratio  is  less  than  10  between 


n  =  0.1  and  0.9.  However,  calculations  of  this  ratio  at  the  same  z  =  0 
location  vary  from  40  (0  <  n  <  0.4)  to  300  (0.4  <  n  <  0.95).  Although 
numerical  diffusion  compounds  the  problem,  the  fact  that  it  is  the  more 
accurate  QUICK  scheme  results  which  show  worse  agreement  with  the  measure 
ments  at  this  location,  strongly  supports  the  contention  that  the  source 
of  false  diffusion  is  more  physical  in  nature  (due  to  the  model)  than  it 
is  numerical  (due  to  the  differencing  scheme).  In  fact,  the  higher 
levels  of  numerical  diffusion  in  the  HYBRID  scheme  distort  physical 
diffusion  in  the  turbulence  model  and  lead  to  the  incorrect  impression 
of  better  predictions  at  various  streamwise  locations! 

In  general,  the  QUICK  scheme  calculations  of  mean  radial  velocity 
show  better  qualitative  agreement  with  the  measurements  than  the  HYBRID 
results.  Both  schemes  yield  calculated  values  of  k  which  are  in  reason¬ 
ably  good  agreement  with  estimates  of  k  from  the  measurements;  with  k  __ 

me  a  S 
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approximated  as  kmeas  =  ^(Ug  +  2  uf).  However,  the  agreement  is  mislead- 

T 

ing.  When  stunned,  the  experimentally  determined  maxima  in  uQ  and  the 

_  u 

2 

minima  in  u£  compensate  to  yield  fairly  uniform  radial  distributions  of 
Sneas'  To  ^ndicate  clearly  the  inadequacy  of  the  model,  the  plots  pro¬ 
vided  show  comparisons  between  predicted  values  of  uQ  and  ur  (calculated 
2  2 

assuming  u  =  u  =  and  direct  measurements  of  these  stresses 

u  r  >3  ca  i  c 

at  2z/Dh  =  0  and  0.5. 

In  several  of  the  plots,  both  of  the  profiles  given  for  U.  at 
2z/Dh  =  0  and  0.5  respectively  show  larger  values  of  this  component  than 
were  actually  measured,  suggesting  different  values  for  the  experimental 
and  computed  mass  flows  in  the  streamwise  direction.  In  fact,  predic¬ 
tions  of  U«  nearer  the  wall,  not  shown  here,  yielded  values  smaller  than 
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those  measured  and  gave  the  required  mass  balance. 

Numerical  experiments  showed  that  a  smaller  length  scale  specifica¬ 
tion  in  equation  (4.4.1)  (i,  =  0.01  Dtl)  at  the  inlet  plane  (Xu  =  -5)  led 

m  H  n 

to  considerably  reduced  levels  of  the  predicted  turbulent  kinetic  energy. 
As  a  consequence,  even  though  there  was  much  worse  agreement  between 
predicted  and  estimated  values  of  k,  in  the  vicinity  of  9  =  90°  profiles 
of  UQ  (but  not  of  11^)  showed  the  local  minima  displayed  by  the  measure¬ 
ments.  Again,  this  type  of  partial  agreement  is  misleading  since  it  is 
artifically  reduced  physical  diffusion  (through  the  specification  of  a 
larger  dissipation  at  the  inlet  plane)  which  yields  the  result.  In 
general,  much  better  overall  predictions  were  obtained  by  prescribing 
an  inlet  plane  length  scale  variations  as  outlined  in  section  4.4 

At  8  3  177°  the  flow  in  the  bend  shows  all  the  maxima  in  the  UQ 
profiles  displaced  towards  the  concave  wall  in  the  presence  of  rela¬ 
tively  large  levels  of  transverse  flow.  Even  though  the  radial  velocity 
component  has  been  decreasing  steadily  from  about  0  =  90°  due  to  the 
oscillatory  inviscid  flow  mechanism,  there  is  no  evidence  in  the  Ur 
profiles  at  any  location  in  the  bend  of  a  Dean  instability  in  the  vicin¬ 
ity  of  the  concave  wall.  The  instability  is  of  the  type  found  by  Taylor 
in  concentric  cylinder  flows  (Dean  [1928])  and  has  been  observed  experi¬ 
mentally  by,  for  example,  Cheng  et  al.  [1979].  It  has  been  predicted 
numerically  only  in  laminar  flow  regime  (Cheng  et  al.  [1975],  Joseph  and 
Smith  [1975],  andGhiaand  Sohkey  [1977]).  Conceivably,  in  a  bend  of 
larger  deflection  angle  than  the  present  one,  the  inviscid  flow  oscilla¬ 
tion  could  also  induce  an  extra  pair  of  vortices  in  the  vicinity  of  the 
concave  wall.  There  is  supporting  evidence  for  this  conjecture  in 


Figure  8  {Station  12)  of  Hawthorne's  paper  [1951]  and  in  the  measured 
results  of  this  work  at  XH  =  1.  At  this  location  the  inviscid  flow 
mechanism  has  not  ceased  to  apply  completely  because  of  the  persistence 
of  streamline  curvature  in  the  flow.  Such  a  mechanism  for  generating  an 
extra  pair  of  vortices  in  the  cross-stream  plane  is  distinct  from  the 
Dean  instability  which  is  associated  with  destabilizing  curvature 
effects.  At  9  =  177°  both  of  the  normal  stresses  display  a  surprising 
degree  of  uniformity  and  the  shear  stress  is  small,  implying  that  the 
flow  is  well-mixed  and  relatively  isotropic  at  this  location.  Except 
for  Ur,  at  this  station  the  calculations  are  in  very  good  agreement 
with  the  measurements.  Although  inaccurate  in  terms  of  absolute  values, 
QUICK  scheme  predictions  of  Ur  are  vastly  superior  to  those  obtained 
with  the  HYBRID  scheme.  The  QUICK  scheme  faithfully  reproduces  the 
measured  profile  curvatures  while  the  HYBRID  scheme  predicts  negative 
values  where  none  exist  in  the  measurements! 

At  XH  =  1  (Figure  6.10)  the  secondary  flow  emerging  from  the  bend 
has  experienced  a  drastic  change  in  both  its  magnitude  and  sense.  The 
radial  velocity  component  has  been  reduced  to  less  than  about  4%  and  a 
transverse  flow  appears  near  the  concave  wall,  opposite  in  sense  to  that 
existing  within  the  bend  at  the  same  radial  location.  A  similar  obser¬ 
vation  was  made  by  Taylor  et  al.  [1982]  in  their  laminar  flow  measure¬ 
ments  at  XH  *  2.5.  As  mentioned  above,  the  phenomenon  is  attributed  to 
the  persistence  of  the  inviscid  flow  mechanism  downstream  of  the  bend. 
Because  of  the  fairly  uniform  and  similar  levels  of  the  normal  stress 
components  measured  at  9  =  177°  and  XH  =  1,  it  is  unlikely  that  gradients 
of  these  stresses  are  altering  significantly  the  intensity  and  pattern 


of  the  cross-stream  flow  at  these  stations. 

As  of  X^|  =  0  the  flow  in  the  downstream  tangent  reverts  slowly  to  a 
straight  duct  flow  condition,  with  the  maximum  in  the  IL  component  tend- 

U 

ing  towards  the  duct  center.  However,  the  maximum  has  moved  only  a 
small  amount,  from  n  =  0.9  to  a  value  of  about  0.65  over  a  duct  length 
of  20  0H-  Similar  indications  of  developing  straight  duct  flow  also 
appear  in  the  profiles  for  the  longitudinal  turbulence  intensity  u0/Ug 
and  the  shear  stress  at  XH  =  5,  10  and  20  (see  Figures  6.11-6.13  and 
6.16).  By  contrast,  the  transverse  components  of  mean  velocity  and  tur¬ 
bulence  intensity  do  not  show  the  same  rates  of  adjustment  to  the  new 
flow  conditions.  The  profiles  for  Ur  indicate  that  as  of  XH  »  5  a 
cross-stream  flow  pattern,  similar  but  much  weaker  to  that  in  the  bend, 
arises  in  the  downstream  tangent.  At  XH  =  2.5,  their  furthest  downstream 
station,  Taylor  et  al.  observed  the  same  result  in  their  turbulent  flow 
measurements. 

The  resurgence  of  a  bend-like  secondary  motion  in  the  downstream 
tangent  is  surprising  and,  at  present,  not  fully  understood.  Calcula¬ 
tions  of  the  flow  between  X^  =  0  and  5,  using  the  QUICK  scheme  are  in 
qualitative  agreement  with  the  measurements.  In  particular,  the  trans¬ 
verse  variations  of  Ur  at  2z/DH  *  0  and  0.5  are  correctly  simulated  by 
QUICK  but  incorrectly  simulated  by  HYBRID.  The  persistence  of  a  peak  in 
UQ  near  the  convex-wall  side  of  the  duct  at  2z/DH  =  0.50  (see  Fig.  6.11) 
is  due  to  turbulence  model  deficiencies.  Because  of  smearing  by  numeri¬ 
cal  diffusion,  this  is  not  reflected  in  the  HYBRID  calculation. 

Predicted  cross-section  vector  plots  and  streamwise  velocity  con¬ 


tours  of  the  flow  at  6  =  177°  using  the  two  differencing  schemes  are 


given  in  Fig.  6.17.  The  results  show  that  the  HYBRID  scheme  predicts 
two  large  vortices  per  half-symmetry  plane  in  the  bend  while  the  QUICK 
scheme  predicts  essentially  one,  with  a  second,  much  smaller  counter¬ 
rotating  vortex,  at  the  convex  wall.  The  HYBRID  results  are  in  blatant 
disagreement  with  the  measurements,  the  QUICK  results  are  in  qualitative 
agreement  and  this  situation  persists  up  to  the  exit  plane  at  =  5. 

6.3  The  ASM  Predictions 

In  this  section  predictions  of  the  flow  using  the  ASM  closure  are 
presented  and  discussed.  The  amount  of  available  core  memory  in  the  CDC 
7600  computer  at  the  Lawrence  Berkeley  Laboratory  required  performing 
the  computations  on  a  grid  different  from  that  reported  in  Section  6.2. 
For  the  present  case,  the  unequally  spaced  grid  distribution  initially 
adopted  was:  (z  =  12)  x  (r  =  22)  x  (x  *  20)  in  the  upstream  tangent  of 
5  Dh;  (z  ■  12)  x  (r  =  22)  x  {9  »  45)  in  the  bend;  and  (z  *  12)  x  (r  =  22) 
x  (x  *  15)  in  the  downstream  tangent  of  4.2  DH.  The  distribution  of 
grid  nodes  over  a  cross-stream  plane  is  given  in  Table  6.2.  Recent  com¬ 
putations  using  a  finer  grid  size  in  the  bend  along  the  main  flow  direc¬ 
tion  yielded  somewhat  better  results.  Accordingly,  60  equally  spaced 
planes  along  the  main  flow  direction  was  distributed  in  the  bend,  a 
1/3  increase  from  what  was  used  in  all  previous  calculations. 

As  mentioned  earlier.  It  has  not  yet  been  possible  to  obtain  stable 
and  converged  results  with  the  QUICK  scheme,  so  that  only  HYBRID  predic¬ 
tions  are  reported  here.  For  comparison  purpose,  new  calculations  using 
the  k-e  model  (HYBRID  only)  with  an  identical  grid  are  also  included. 
Transverse  distributions  of  the  predictions  at  z/DH  «  0  and  z/DH  *  0.25 


for  several  streamwise  planes  are  plotted  in  Figures  6.18  to  6.37.  In 
these  figures,  solid  lines  represent  best  curve  fits  to  the  experimental 
data  discussed  earlier,  dot-dash  lines  are  the  ASM  predictions,  and  dot- 
dot-dash  lines  are  the  corresponding  k-e  model  predictions. 

From  an  overall  inspection  of  the  results,  several  observations 
concerning  the  ASM  predictions  can  be  made.  They  are: 

1.  The  mean  flow  results  are  not  better  than  those  predicted  with 
the  k-e  model.  In  fact,  at  9  =  90°  (Figure  6.22),  it  is  the  k-e  model 
which  produces  better  agreement  with  the  data! 

2.  The  ASM  predicts  a  slower  development  of  the  secondary  motion 
in  the  bend  than  that  predicted  by  the  k-e  model. 

3.  Although  capable  of  predicting  turbulence  anistropy  as  shown 

in  the  plots  of  ur/Ug  and  Q0/Ug,  the  ASM  predictions  failed  to  reproduce 

the  experimentally  observed  local  extrema  in  Ua,  U  ,  uQ,  0  and  u  uQ 

u  r  u  r  r  u 

around  n  =  0.4  and  0  =  90°. 

Since  the  equation  system  describing  the  flow  is  highly  non-linear 
and  strongly  coupled,  it  is  very  difficult  to  single  out  the  causes  for 
the  poor  agreement  revealed  in  Figures  6.18  to  6.37.  Faced  with  these 
discouraging  results,  one  might  well  ask  whether  or  not  the  ASM  closure 
proposed  in  Section  3.3.2,  which  was  developed  with  2D  flows  in  mind, 
is  capable  of  predicting  such  a  complex,  three-dimensional,  highly 
anisotropic  turbulent  flow  as  the  present  one.  To  answer  this  question, 
or  to  at  least  bear  light  on  the  problem,  two  simplified  "experimental" 
tests  were  conducted  as  explained  below. 

First,  the  measurements  of  Up,  U0  along  the  symmetry  plane  at 
9  *  90°  (where  the  predictions  were  the  worst),  together  with  the  k-c 


model  estimates  of  k,  e,  and  p  were  substituted  into  the  ASM  relations 

given  in  Section  3.3.3.*  Solving  the  resultant  equation  system  yielded 

the  distributions  plotted  in  Figures  6.23,  6.25  and  6.38  as  the  dashed 

lines.  It  is  seen  that,  except  for  ur,  the  ASM  is  capable  of  predicting 

the  correct  trends  of  the  flow  anistropy  provided  that  accurate  spatial 

resolution  of  the  mean  velocities  is  presented  to  the  equation  system. 

A  second,  more  direct  test  was  also  performed.  In  addition  to  Ur 

and  U. ,  the  measurements  of  u_  and  CL  were  also  used.  Substituting  these 
y  r  y 

measured  quantities,  together  with  the  k-e  model  estimates  of  k,  e,  and 

P  into  equation  (3. 3. 3. 4)  and  solving  for  uTuT  yielded  the  much  better 

r  y 

distribution  plotted  in  Figure  6.38  as  the  dotted  line. 

As  a  result  of  these  two  tests,  it  is  now  believed  that  the  ASM 
closure  proposed  in  Section  3.3.2  is  indeed  adequate  for  predicting  com¬ 
plex  anistropic- 3D  turbulent  flows  of  the  present  kind,  provided  that 
accurate  spatial  resolution  of  the  mean  velocities  is  achieved.  Previ ous 
calculations  using  the  k-e  model  suggested,  however,  that  substantially 
refined  grid  would  be  required  to  significantly  improve  the  numerical 
accuracy  of  the  present  predictions.  Computation  costs  and  memory  avail¬ 
ability  simply  prohibited  the  use  of  finer  grids  to  reduce  the  numerical 
diffusion.  The  alternative  of  using  higher  order  differencing  such  as 
the  QUICK  scheme  is  probably  a  better  course  to  pursue. 

During  the  course  of  the  research,  three  variations  of  the  log-wall 
function  approach  were  tested:  the  standard  TEACH  code  method  (Gosman 

*In  the  absence  of  information  concerning  streamwise  and  spanwlse  varia¬ 
tions  of  the  velocity  components,  terms  containing  these  quantities, 

including  3UQ/30,  3U.730  and  3U„/3z,  were  neglected, 
y  r  z 


and  Ideriah  [1976]);  the  Chieng  and  Launder  [1979]  method;  and  the 
method  used  by,  among  others,  Naot  and  Rodi  [1981].  As  a  consequence, 
it  was  found  that  the  flow  pattern,  especially  the  secondary  flow 
velocity  components  Ur  and  U2,  is  sensitive  to  the  way  the  near  wall 
flow  is  treated.  However,  no  one  of  the  three  methods  produced  results 
in  better  overal  agreement  with  the  data  than  the  others.  Because  of 
its  simplicity,  the  treatment  of  Naot  and  Rodi  was  chosen  for  the  cal¬ 
culations  presented  above. 


CHAPTER  7.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  new  measurements  provided  here  for  the  mean  flow  in  strongly 
curved  180°  bends  reveal  features  which  are  in  qualitative  agreement 
with  results  obtained  from  inviscid  flow  analysis.  Measurements  of  the 
turbulence  characteristics  of  the  flow  at  9  =  90°  and  130°  show  striking 
variations  in  the  radial  direction  which  appear  to  be  due  to  large 
shearing  motions  induced  by  the  inviscid  mechanism  in  the  core  of  the 
flow.  Within  the  context  of  an  interpretation  which  neglects  turbulence 
diffusion  and  pressure  redistribution  effects,  the  velocity  gradients 
induced  in  the  core  of  the  flow  work  locally  on  the  turbulent  stresses 
to  increase  the  longitudinal  stress  component  (uQ)  and  decrease  the 
radial  stress  component  (uf). 

Between  e  =  0°  and  45°,  stabilizing  and  destabilizing  curvature 
effects  at  the  convex  and  concave  wall  respectively  are  responsible  for 
generating  relatively  large  levels  of  anisotropy  near  these  surfaces. 

As  of  9  «  45°,  the  cross-stream  motion  contributes  to  the  production  of 
a  more  complicated  anisotropic  pattern  in  the  turbulent  stresses.  For 
example,  regions  of  negative  production  of  turbulent  kinetic  energy  arise 
at  0  *  90°.  Between  9  *  130°  and  177°  the  anistropy  in  the  measured 
normal  stresses  is  reduced,  due  to  the  uni formi zing  influence  of  the 
secondary  motion  in  the  flow. 

Measurements  of  Ur  In  the  bend  show  that,  except  for  around  6  *  45°, 
the  return  flow  in  the  cross-stream  plane  is  restricted  to  fairly  narrow 
bands  of  width  »  DH/8  adjacent  to  each  of  the  side  walls.  The  measure- 
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ments  of  U  did  not  reveal  vortex  structures  at  the  convex  wall  in  the 
r 

bend,  nor  at  the  concave  wall  where  they  might  be  expected  to  arise  due 
to  the  Dean  instability.  However,  the  sense  of  the  cross-stream  flow 
at  XH  =  1  does  suggest  a  brief  appearance  of  this  type  of  instability  in 
the  exit  tangent.  Between  6=0°  and  X^  =  20  there  was  no  evidence  in 
the  measurements  of  a  significant  Reynolds  stress  driven  secondary 
motion. 

In  this  investigation  the  influence  of  the  bend  on  the  upstream 
tangent  flow  is  already  noticeable  at  XH  =  -5.  Although  weak,  a  trans¬ 
verse  flow  is  set  up  at  this  location  (and  at  X^  =  -1)  of  characteris¬ 
tics  which  correspond  with  the  simple  notion  of  a  superposition  of  two 
cross-stream  motions:  the  first  induced  by  Reynolds  stress  gradients 
in  the  cross-stream  plane;  the  second  due  to  the  favorable  streamwise 
pressure  gradient  arising  at  the  convex  wall  in  the  bend.  In  the  down¬ 
stream  tangent  secondary  motions  are  drastically  decreased  between  e  = 

180°  and  XH  =  1 .  The  decrease  is  related  to  the  oscillating  nature  of 
the  secondary  motion  which,  as  shown  by  inviscid  flow  theory  and  supported 
by  the  present  results,  as  of  9  =  90°  has  started  opposing  the  initial 
source  of  circulation  set  up  by  the  transverse  pressure  gradient  in  the 
bend.  Between  X^  ■  1  and  20  the  flow  reestablishes  characteristics 
typical  of  straight  ducts.  However,  the  measurements  show  the  persis¬ 
tence  of  a  weak  secondary  motion  in  the  downstream  tangent  with  the  same 
sense  of  rotation  as  the  flow  in  the  bend  and  a  slow  recovery  of  the 
radial  normal  stress  compared  to  the  longitudinal  component. 

Predictions  of  the  experimental  flow  were  first  made  with  a  k-e 
model  of  turbulence  using  a  partially-parabolic  numerical  procedure.  Two 


finite  differencing  techniques  were  used  with  global  accuracies  of  first 
order  (HYBRID)  and  second  order  (QUICK).  Testing  and  predictions  with 
QUICK  in  the  laminar  flow  regime  clearly  showed  it  to  be  superior  to 
HYBRID.  The  differences  between  the  results  obtained  using  the  two 
schemes  in  turbulent  flow  are  notable.  Of  the  two,  only  QUICK  repro¬ 
duces  correctly  the  trends  in  the  measured  cross-stream  flow.  The 
HYBRID  scheme  gives  the  incorrect  impression  of  better  streamwise  com¬ 
ponent  velocity  predictions.  This  is  due  to  numerical  diffusion  which 
smears  out  the  spacial  variation  of  this  component.  Although  QUICK 
scheme  predictions  for  Ua  differ  more  from  the  measurements,  this  is 
explained  by  the  fact  that  with  numerical  diffusion  reduced,  it  is  the 
false  physical  diffusion  in  the  k-e  model  which  is  revealed  by  the  com¬ 
puted  results. 

Numerical  experiments  confirmed  a  fairly  strong  sensitivity  of  the 
model  to  the  inlet  plane  specification  of  dissipation.  Fixing  a  small- 
but  incorrect  dissipation  length  scale  at  the  inlet  plane  led  to  some¬ 
what  improved  predictions  of  the  UQ  component  at  0  =  90°,  where  this 
component  shows  large  changes  in  the  transverse  direction.  However, 
predictions  of  Uf  and  of  k  were  much  less  satisfactory. 

In  an  earlier  study  (Chang  et  al.  [1982])  it  was  reported  that 
predictions  of  turbulent  flow  in  a  90°  bend  were  mildly  sensitive  to 
the  sense  of  cross-stream  flow  fixed  at  the  inlet  boundary  condition. 

The  effect  was  too  small  to  alter  significantly  the  streamwise  evolution 
of  the  flow.  In  this  study,  insignificant  differences  were  observed 
between  predictions  of  the  flow  with  the  cross-stream  motion  set  equal 
to  zero,  and  predictions  with  the  cross-stream  motion  specified  from 


the  measurements.  The  insensitivity  of  the  flow  in  a  bend  to  the  inlet 
plane  streamwise  component  of  vorticity  is  explained  by  reference  to 
the  general  form  of  the  equation  for  the  variation  of  this  quantity 
along  a  streamline  (Horlock  and  Lakshminarayana  [1973]): 


In  the  case  of  a  flow  entering  a  bend  from  a  straight  duct,  generally 

£  >>  >  especially  near  the  side  walls.  As  a  consequence,  the  equa- 

r  u 

tion  shows  that  along  streamlines  curving  through  the  bend,  the  stream- 
wise  development  of  vorticity  is  due  primarily  to  deflection  of  trans¬ 
verse  vorticity.  Viscous  diffusion  of  streamwise  vorticity,  even  in 
turbulent  flows,  is  of  secondary  importance. 

Subsequent  predictions  with  the  ASM  closure  did  not  yield  better 
agreement  with  measurements  than  did  the  k-e  results.  Consequently  , 
two  "experimental"  tests  were  carried  out  to  check  indirectly  the  pre¬ 
dicting  ability  of  the  ASM  closure.  It  was  found  that  the  approximations 
proposed  in  Section  3.3.2  are  indeed  capable  of  resolving  the  anisotropy 
of  the  turbulence  correctly,  provided  that  the  mean  velocity  field  is 
known  accurately;  i.e.  with  good  spatial  resolution. 

The  following  recommendations  are  offered  for  continuing  work: 

1.  Laminar  flow  measurements  should  be  performed  in  the  bend  and 
downstream  tangent  to  establish  definitively  the  superior  perfor¬ 
mance  of  QUICK  for  predicting  curved  duct  flows  and  for  verifying 
the  various  vortical  structures  predicted  by  this  scheme. 

2.  More  detailed  turbulent  flow  measurements  are  required  to  under- 
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stand  the  relaxation  processes  taking  place  in  the  downstream  tan¬ 
gent. 

3.  Modeling  at  the  full  stress  equation  closure  level  will  certainly 
be  required  to  predict  wall  interactions  accurately,  including 
curvature  effects. 

4.  Fine  grid  calculation  with  the  QUICK  scheme  using  the  ASM  closure 
is  recommended.  This  could  allow  more  definite  conclusions  to  be 
drawn  concerning  the  behavior  of  the  ASM  closure. 

5.  It  is  believed  that  the  present  specifications  of  the  log-wall 
relations  along  all  three  walls  are  inaccurate.  Since  they  are 
essential  to  produce  reliable  predictions  by  means  of  schemes 
which  invoke  these  relations,  more  detailed  experimental  work 
concerning  these  laws  in  3D  turbulent  walls  flows  is  strongly 
recommended. 


6.  In  parallel  to  the  experimental  work  mentioned  above  in  5., 

numerical  experiments  on  the  sensitivity  of  the  flow  to  the  model 
constants  <  and  E  are  recommended.  In  particular,  the  constant  E 
should  be  made  functions  of  pressure  gradients  along  the  main  flow 
direction  for  Ufl  and  over  a  cross-sectional  plane  for  U„.  and  U  . 
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Quantity 

Range  of  Maximum 
Systematic  Error  + 
(deviation  5) 

Maximum 
Random  Error 
(r.m.s.  5) 

VUB 

-  15  to  25 

±  25 

VUB 

-  15  to  25 

±  25  to  ±  35 

VUB 

-  25  to  25 

±  25  to  ±  35 

VUB 

-  25  to  25 

±  25  to  ±  45 

VVUB 

-  25  to  25 

±  35  to  ±  75 

+  Error  confined  mainly  to  near  wall  flow  regions; 

does  not  Include  angular  uncertainty  discussed 
In  text. 

++  Largest  errors  confined  to  small  flow  regions;  does 

not  Include  a  small  positioning  uncertainty  of  probe 
volume  In  the  flow. 


Table  2.1  Estimated  maximum  measurement  errors 
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T?  i»  6.1  -rid  distribution  over  a  cross-stream  plane  for  the  k-e 
model  calculations 
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Figure  1.1  General  flow  pattern  of  the  secondary  motion  of  the  first 
kind  in  a  curved  duct  of  square  cross-section. 


ure  2.2  Optimum  experimental  arrangement  of  stainless  steel  screens 
and  perforated  plate  In  the  flow-unlformizlng  section 
preceding  the  upstream  tangent. 


Figure  2.3  Top  and  side  view  of  laser-Doppler  velocimeter  aligned  for 
measurements  at  a  bend  angle  6  «  180°.  Traversing  mechanism 
details  are  not  shown. 
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Fully  developed  laminar  flow  in  a  straight 
wall  moving  at  a  constant  speed  normal  to 
direction. 


Figure  5.5  k-e  model  calculation  of  fully-developed  turbulent  flow  i 
a  20  straight  channel;  mean  longitudinal  velocity  U/U 


Figure  5.7  Fully  developed  turbulent  flow  In  a  2D  straight  channel; 

distributions  of  u  ,  u  and  u  calculated  by  using  the  ASM 
relations  and  U  ,  x  k  “and  e  fiom  the  k-e  model  calculation 
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Figure  5.8  ASM  calculation  of  fully  developed  turbulent  flow  in 

straight  channel;  mean  longitudinal  velocity  U/U  . 

max 


(QS60  vi# q  s#Minvn 


Figure  5.11  ASM  calculation  of  fully  developed  turbulent  flow  In  a  2D 
curved  channel;  uQ/Ueand  U0/Ue  max:  ( — )  measurements  of 

Esklnazl  and  Yeh  [1956];  (-•-)  ASM  calculation;  ( - )  ASM 

calculation  with  non-symmetrlc  f-functlon  of  Humphrey  and 
Pourahmadl  [1982]. 


Figure  5.12  ASM  calculation  of  fully  developed  turbulent  flow  in  a  2D 

curved  channel;  uua/Un  and  tT/Ua:  ( —  and  o)  measure- 

r  w  o  >max  r  o 

ments;  (-•-)  ASM  calculation;  ASM  calculation  with 

non-symmetri c  f-function  of  Humphrey  and  Pourahmadi  [1982]. 


Figure  6.1  Distribution  of  pressure  coefficient  (Cp  *  AP/pUg)  in  the 

bend  and  tangents  for  conditions  of  the  experiment:  Re  * 
56,700,  Rc/Dh  *  3.35  and  De  *  21,900.  Measurements  made 

through  pressure  taps  on  side  walls  of  test  section  compon 
s. 


05 

r»* 


•#OrtO®o9#®00® 


*0  _  O9o0 


O««0o° 

. 


Oo°oo 

oo®»°  eo« 


0*  10 


VU1.4.1100 


0  OS 


(Sr/lWOO 


OS 


Figure  6.2  Comparison  between  turbulent  flow  measurements  of  mean 

velocity  and  normal  stress  from  this  work  (o)  at  XH  *  -1  a 

measurements  by  Mel  ling  and  White! aw  ( - )  after  a  devel 

opment  length  of  36.8  Du. 
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6.3  Turbulent  flow  measurements  of  mean  velocity  and  normal 
stress  in  the  upstream  tangent  at  XH  »  -5. 
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Figure  6.4  Turbulent  measurements  of  mean  velocity  and  normal  stress 
In  the  upstream  tangent  at  Xu  *  -1. 
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Figure  6.7  Turbulent  flow  measurements  and  k-e  model  calculations  of 
mean  velocity  and  normal  stress  In  the  bend  at  e  *  90°: 

( - )  HYBRID,  ( — )  QUICK;  dashed  line  at  0.875  is  best 

fit  to  measurements. 


Turbulent  flow  measurements  and  k-e  model  calculati 
mean  velocity  and  normal  stress  In  the  bend  at  6  = 
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Figure  6.9  Turbulent  flow  measurements  and  k-e  model  calculations  of 
mean  velocity  and  normal  stress  In  the  bend  at  9  *  177°; 

( - )  HYBRID,  ( — )  QUICK;  dashed  lines  at  0.75  and  0.875 

are  best  fits  to  measurements. 
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Figure  6.10  Turbulent  flow  measurements  of  mean  velocity  and  normal 
stress  in  the  downstream  tangent  at  Xu  «  1. 
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Figure  6.11  Turbulent  flow  measurements  and  k-e  model  calculations  of 
mean  velocity  and  normal  stress  In  the  downstream  tangent 
at  XH  «  5:  ( - )  HYBRID,  ( — )  QUICK;  dashed  line  at  0.875 

Is  best  to  fit  measurements. 
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Figure  6.12  Turbulent  flow  measurements  of  mean  velocity  and  normal 
stress  In  the  downstream  tangent  at  XH  «  10;  dashed  lines 

at  0.75  and  0.875  are  best  fits  to  measurements. 
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Figure  6.13  Turbulent  flow  measurements  of  mean  velocity  and  normal 
stress  In  the  downstream  tangent  at  XH  ■  20;  dashed  line 

at  0.875  is  best  fit  to  measurements. 
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Figure  6.15  Turbulent  flow  measurements  of  shear  stress  at  four  longi 
tudlnal  stations  In  the  bend;  dashed  line  at  0.75  Is  best 
fit  to  measurements. 


6.19  ASM  calculations  of  (Qg/Ug)  x  100  at  6 
fits  to  measurements,  (-•-)  ASM. 
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Figure  6.22  ASM  calculations  of  U6/Ug  at  6  •  90°:  ( — )  best  fits 

to  measurements,  (-•-)  ASM,  (-••-)  k-e  model. 


Figure  6.25  ASM  calculations  of  (ur/Ug)  x  100  at  6  *  90°:  ( — )  best 

fits  to  measurements,  (-•-)  ASM,  ( — )  "experimental" 
test  #1. 
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Figure  6.28  ASM  calculations  of  (Ur/Ug)  at  e  •  177#:  (  )  best  fits 

to  measurements,  (-•-)  ASM,  (-•— )  k-e  model. 
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Figure  6.31  ASM  calculations  of  (uQ/Ug)  x  100  at  x/D 
fits  to  measurements,  (-•-)  ASM. 


